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The Oscillation and Noise of an Overpressure 
Sonic Jet! 


ANDREW G. HAMMITT* 
General Applied Science Laboratories, Inc. 


Summary 


A series of experiments on plane, two-dimensional, overpressure 
sonic jets has been performed to study the oscillating behavior of 
the jet and accompanying sound field. Shadowgraphs were ob- 
tained which show the oscillations and generation of sound waves 
in some detail. Measurements from these pictures show that 
the oscillating jet and sound fields can be described in terms of 
the length of the shock-wave cells. The importance of the sound 
waves acting upon the base of the jet is demonstrated and it is 
found that the jet can be stabilized by shielding its base from 


these sound waves. 


Symbols 
speed of sound in atmosphere around jet 
constants of proportionality 


width of jet 

distance from origin of jet to source of sound waves 
length of shock-wave cells 

Mach number 

frequency of oscillation 

atmospheric pressure 

stagnation pressure to cause critical flow 
stagnation pressure 

Strouhal number 

wavelength of jet oscillations 

wavelength of sound waves 


Introduction 


A HIGH-VELOCITY JET OF AIR has been observed to 
be a source of considerable noise. In general, 
this noise covers a broad range of frequencies. For 

Received May 11, 1960. Revised and received August 2, 
1960. 
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the particular case of an overpressure jet issuing from 
a sonic orifice, a single-frequency tone is sometimes 
heard which predominates over the other noise. This 
generation of a pure tone has been called jet screech. 
The purpose of this study is to obtain fundamental 
information on the mechanism of the jet screech phe- 


nomenon. 


Experimental Investigation 


The experimental investigation to be described was 
carried out at Princeton University using the facilities 
of the Gas Dynamics Laboratory. A settling chamber 
was built to provide a low-turbulence, low-noise-level 
supply for various jet orifices during operation from 
an available 3,000-psi air supply. Since maximum 
pressures in the settling chamber of the order of 100 psi 
were desired, a mechanism had to be provided to 
throttle the 3,000-psi air supply without the introduc- 
tion of turbulence or noise. To accomplish this throt- 
tling, a valve followed by two pieces of porous stainless 
steel was used. A drawing of this settling chamber 
together with the nozzles used is shown in Fig. 1. Two 
plane, two-dimensional slit-type jet nozzles were used, 
one with the jet originating at the face of the flange 
and the other with a tube or spout so that the jet 
originates away from the face of the flange. Both 
nozzles were of the converging type. 

Studies of the jets were made using flow visualization 
techniques and microphone pressure measurements 
of the sound field. The shadowgraph technique for 
visualizing the second derivative of the density field was 
used extensively. The sensitivity of this system can 
be adjusted by placing the film at different distances 
from the phenomena being observed, the greater the 
distance the more sensitive is the system. Pictures 
near the jet show the stronger phenomena more clearly 
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Fic. l(a). Settling chamber used in jet-noise studies. 


while those further away emphasize the weaker phe- 
nomena. A high-voltage spark with a duration of 
about 1 microsecond was used as the light source, and 
proved adequately fast to stop the motion. 

The microphones used were a Massa Type M-141 
crystal microphone having a frequency range to 30 ke 
and a General Radio Rochelle-salts microphone with a 
frequency range to 10 ke. The Massa microphone 
was used for obtaining all data presented except that 
shown in Fig. 8. In this latter case, the higher response 
of the Rochelle-salts microphone was required. 


Experimental Results 


Optical Studies 


A series of shadowgraphs of the nozzle with spout 
taken looking at the narrow edge of the jet, is shown in 
Fig. 2. These shadowgraphs were taken with the film 
at two distances from the jet, 5 in. and 4 ft. The 5-in. 
pictures show the jet structure, while the 4-ft pictures 
show the sound field. The shadowgraphs shown are 
for a range of stagnation pressures varying from 30 to 
64 psia. 

At 30 psia just the 5-in. picture is shown, since the 
jet is stable at this condition. In the object just to 
the right of the jet are two small holes 1 in. apart 
which furnish a reference scale for these pictures. The 
object appearing on the right side of the picture is a 
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microphone. Fig. 1(b) shows the jet at 39 psia from 
both the near and the far positions. At this pressure, 
the jet is violently unstable. The closeup shows the 
oscillations of the jet structure, and the distant view 
shows the sound waves spreading away from the jet. 
It is interesting to note the much more rapid spread 
of the jet at this pressure than at the 30 psia condition. 
This rapid spreading is characteristic of the oscillating 
jet. Figs. l(c) at 54 psia and 1(d) at 64 psia show the 
same phenomena as 1(b) except that the sound waves 
are stronger and of lower frequency. Figs. l(c) and 
l(b) also show a series of waves with small curvature 
lying in the horizontal plane. These waves are the re- 
flections from the flange on the settling chamber at 
the base of the nozzle. 

The predominant waves originating from the jet 
can be seen to be approximately semicircular and to 
originate alternately from each side of the jet. The 
sources of these waves on each side of the jet seem to be 
located at about the same height. A detailed ex- 
amination of the pictures presented here, and of others 
as well shows that apparently some weaker waves 
originate from other heights also, in the same alter- 
nating fashion. An examination of the pictures of the 
jet structure shows that the regular jet shock pattern 
can be seen through the third cell and that the jet is rela- 
tively stable up to this point. After the third cell, 
the violent oscillations take place and the predominant 
sound waves are emitted near the top of the third cell. 

Pictures of this oscillating jet at various times during 
a complete cycle would be useful in determining the 
manner in which the sound waves originate. The 
jet is oscillating at between 5,000 and 10,000 cycles 
per second, and since several pictures per cycle are 
necessary, a frame speed of about 50,000 frames per 
second would be required. A camera of this speed 
was not available, and so a different technique was 
used. A series of pictures was taken of the jet at one 
stagnation pressure. These pictures, having their own 
built-in clock in the position of the sound waves, could 
then be arranged in chronological order. Fig. 3 shows 
a series of such pictures taken at Py = 64 psia and 
arranged in this manner. The time intervals between 
the pictures can be judged by the progress of the sound 
waves. The stability of the oscillating process is demon- 
strated by the fact that these pictures, taken millions 
of cycles apart, can be considered to represent a single 
cycle. The pictures show that the oscillations, which 
are small near the base of the jet, build up as they 


Fic. 1(b) Left: Plane two-dimensional nozzle with spout. 1(c) 
Right: Plane two-dimensional nozzle without spout. 
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Fic.2. Shadowgraphs of jet at various P». 





Top: film located 5in. from jet; bottom: film located 4 ft. from jet 


(a) Py = 30 psia; (b) Po = 39 psia; (c) Po 54 psia; (d) Po 64 psia 


travel out along the jet. After the third cell vortex 
tails appear which are left behind in an alternating 
manner. A new sound wave seems to originate when 
the jet returns to and contacts this vortex tail. These 
observations are interesting but offer no real explana- 
tion of the origin of the sound waves at this point. 

Some quantitative observations also can be made 
from these pictures. The size of the shock-wave 
cells was measured and the results from these and 
similar pictures are shown in Fig. 4 as a function of 
stagnation pressure. The wavelength of the sound 
waves, the wavelengths of the jet oscillations, and the 
location of the apparent source of the sound waves can 
all be determined. At some pressures, it was possible 
to locate the source of both the major and minor wave 
system, and both of these are shown in Fig. 4. These 
results compare favorably with similar measurements 
on both the cell size and sound wavelength.' It should 
be pointed out now, and it will be shown in greater de- 
tail later, that the oscillations do not grow uniformly 
stronger with stagnation pressure but reach a peak, 
fall off, and then rise to a new peak. Between these 
peak regions, the frequency can be somewhat confused. 
The wavelengths may go in steps with increasing stag- 
nation pressures, as indicated in Fig. 4. The pictures 
chosen to be presented are those at the peaks of the 
oscillations, since they are the most interesting. 

If the sound wavelength, the jet oscillation wave- 
length, and the distance to the source of the sound are 


all nondimensionalized by the shock-cell size, the re- 
sults shown in Fig. 5 are found. It would appear 
that all these quantities have the same dependence on 
the stagnation pressure, and all are related to each other 


by constants. 


Micropnone Studies 

A more complete analysis of the frequency and 
amplitude of the sound waves as a function of stagna- 
tion pressure was obtained by analyzing the output 
of the microphone when located beside the jet in the 
position shown in Figs. 2(a) and 3(a). The results 
are shown in Figs. 6(a) and 6(b). The frequency 
decreases continuously as the stagnation pressure in- 
creases, and the amplitude has several peaks which 
suggest that the curve be broken into separate regions. 
These experiments were made in two different ways. 
As observed in the shadowgraphs (Fig. 2), there are 
sound waves reflected from the base flange entering 
the sound field around the jet. One set of measure- 
ments was taken with this configuration, and one with 
a sound-deadening material placed on the base flange 
so that the reflected waves were eliminated. The 
general frequency—stagnation pressure curve is sub- 
stantially unaffected by this change but there are 
definite changes in the locations of the peak ampli- 
tudes. The frequency jumps that seem to take place 
between the regions of peak amplitude are decreased 
and almost eliminated with the use of the sound-absorb- 
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Fic. 3(a). Shadowgraphs of jet at 64 psia. 
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ing material. The frequency measured without sound- 
absorbing material on the flange sometimes showed an 
overlap when the frequency jump occurred. It is 
believed that this effect is caused by the oscillation 
being unstable and switching between the two modes. 
Pictures of the jet taken in these regions suggest this 
explanation, but complete proof is not available.  Pic- 
tures of the jet both with and without the sound- 
absorbing material on the flange show the same results 
as the microphone records. These results suggest 
that the jet oscillations“are sensitive to the external 
sound waves in the vicinity of the jet. 





5 


Fic. 3(b). Shadowgraphs of jet at 64 psia. 


Series show chronological history of jet oscillations. 


The direction of the sound radiation was also studied 
by means of the microphone. The shadowgraphs of 
the jet show a very prominent set of sound waves 
radiating from the jet in a circular shape. The micro- 
phone was mounted in a quadrant with the center 
located at the source of this sound-wave system, and 
measurements of wave forms were taken at various 
angular positions. Photographs of the wave forms 
as shown on an oscillograph are presented in Fig. 7. 
The angular positions are measured from the jet axis 
with 0° indicating the upstream direction. As the 
microphone approaches the downstream direction (140 
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OVERPRESSURE SONIC JET 


and 160°) the picture becomes more fuzzy, indicating 
the increase in background noise field and possible direct 
jet impingement. The 10° and 160° positions indj- 
cate an approximately sinusoidal variation at the 
fundamental frequency, corresponding to the wave- 
lengths of the principal wave system. The 90° posi- 
tion shows another sinusoidal wave picture but at 
double the frequency of the fundamental system. 
The other pictures show the intermediate steps in 
changing from one system to the other. 

These results, while not suggested by the shadow- 
graphs, are consistent with results presented in reference 
1. The explanation suggested in reference | is that 
the interference of the sound waves from a series of 
sources located along the jet causes this effect. If the 
additional source, which is visible in some of the 
shadowgraphs, is considered, it can be determined that 
its location and phase are such that sound waves from 
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the two sources reinforce each other in the 0° and 
180° directions and interfere in the 90° direction. In 
the 90° direction the second harmonic of each source 
will reinforce and cause the double-frequency response 
observed. The amplitude shown at the 90° position 
gives a measure of the combined second harmonic com- 
ponent of these waves. The sound waves from these 
various sources would have to be of about the same 
strength to cause this complete interference of the first 
harmonic. This conclusion is inconsistent with the 
relative strength shown on the shadowgraphs. 


Effect of Sound Waves Impinging on the Base of the Jet 


References 1 and 2 suggest that the sound waves in 
the external field act upon the jet near its base to cause 
disturbances which are amplified by the jet and which, 
in turn, cause new sound waves. The sensitivity of 
low-speed laminar jets to sound waves has been well 
The mechanism suggested in references 
The effects observed 


established.’ 
| and 2 is a feedback mechanism. 
in these experiments due to change of the sound re- 
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flection condition on the base flange lend support to 
these ideas. It follows from this reasoning that the 
jet could be stabilized if this feedback loop were 
broken. 

The method that was used to break the loop was a 
sound-deadening and reflecting baifle alongside the jet 
between the source of the sound and the jet base (Fig. 
8). This mechanism provides a means of adjusting the 
vertical position of the two baifle plates and the size 
of the slit between them. Experiments were run in 
which the baffle plates served as reflecting surfaces 
and, when covered with soft material, as absorbing 
surfaces. 
the sound waves from propagating from their source to 
the base of the jet but not to interfere with the jet-flow 
field. While these two requirements are incompatible, 
it does appear possible to effectively block the return- 


The purpose of these baffles is to prevent 


ing sound waves without extending the baffles into the 
high-velocity part of the jet. There is some interfer- 
ence with the external flow field induced by the jet, 
but this interference is probably not very important. 
By use of these baffles, it was possible to stabilize the 
jet oscillations. Typical pictures with and without 
the baffles are shown in Fig. 9. 
jet oscillations without the baifles. 
tions take place beginning at the end of the third cell. 
In Fig. 9(b) the baffles are used with a piece of sound- 


Fig. 9(a) shows the 
Definite oscilla- 
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absorbing material on base flange as a function of Py 


absorbing material about '/2 in. thick and | in. high on 
each side of the jet, the metal on the top surface of the 
baffle being bare. The shadow of the fringes of the 
sound-absorbing material as blown up by the jet can 
be seen in the picture. With this baffle arrangement, 
the jet is stabilized so that there is almost no detectable 
oscillations, even over an extended field of vision. 
The two figures are presented so that the jet nozzle is 
in the same position in each picture. 

Several configurations with the baffles were tried with 
different use of sound-absorbing material. The jet 
was stabilized by several of these arrangements. The 
simplest arrangement was obtained by using the flat 
slit nozzle and attaching sound-absorbing material on 
the face. The sound-absorbing material was cut back 
so that it could interfere very little with the jet flow. 
This configuration was very successful in stabilizing 
the jet and eliminating the characteristic screech. The 
flow with this configuration, both with and without the 
absorbing material, is shown in Figs. 10(a) and 10(b). 
This small change had an amazing effect. This 
arrangement appears to be an effective practical 


ane 


Fic. 7. Pictures of oscilloscope trace for microphone placed at 
different positions rotating about source of principal sound waves. 
P) = 64 psia. Angles measured from jet-axis upsteam direction: 
(a) 10°; (b) 60°; (c) 90°; (d) 120°; (e) 140°; (f) 160°. 
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means for prevention of jet oscillations of this type. 

While the stabilization of the jet caused almost 
complete elimination of the characteristic screech fre- 
quency, it does not affect the background roar of the 
jet. The measurements of the overall jet noise show 
that it is reduced by only about 3 db out of 100 db, how- 
ever, the tone changes from the characteristic screech 
to a roar. 

Since the feedback loop had been successfully broken 
by the use of the baffles, an attempt was made to pro- 
vide a substitute means of excitation. For this pur- 
pose, a 40-watt speaker was located so as to cause 
sound waves to impinge on the base of the jet (Fig. 8). 
The speaker could not exactly reproduce the sound 
waves cut off by the baffles, but it was the best arrange- 
ment that could be made simply. Tests were run 
with the speaker driven at a given frequency and a 
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Fic. 8. Jet-baffle system with speaker 


microphone located above the baifle. The output of 
the microphone was connected to a frequency analyzer 
tuned to the oscillator frequency. The output of the 
analyzer was then recorded as a function of jet stagna- 
The results at three different oscillator 
At about 42 psia 


tion pressure. 
frequencies are shown in Fig. 11. 
stagnation pressure the microphone output peaks at 
all frequencies, but there is little effect of the speaker. 
At higher pressures, there is another noise peak that 
is very much affected by the speaker output. The 
stagnation pressures and frequencies at which this peak 
occurs are roughly those shown in Table 1. The last 
column shows the stagnation pressures at which nat- 
ural vibrations at the three frequencies were ob- 
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Shadowgraph of two-dimensional nozzle at Py = 48 psia: 


(a) without baffles; (b) with baffles 


served, taken from Fig. 6. There is a fair correspond- 
ence between the last two columns, showing that the 
jet amplifies the sound at particular frequencies corre- 
sponding roughly to the frequencies observed for the 


jet screech. 


Discussion 


In references | and 2, the suggestion is made that the 
frequency-determining mechanism may involve the 
tuning of the feedback system. It is then shown that 
(l/X) + (//n) = const. + integer, and this result is 
cited as a proof of the feedback mechanism. The pres- 
ent results are in agreement with this timing mecha- 
nism. In Fig. 5, it was shown that all the quantities 
l, \, 7 were related to L, the shock-cell size, by con- 
stants. The quantity /, which in some pictures was 
found to have multiple values, seems to be given by 
the formula //L = const. + (0 or 1). From this rela- 
tion to the shock-cell size, it is obvious that (//A) + 
(1/y) = const. for either value of //L. The fact that 
this sum differs by an integer when // Z changes by | is 
the result of the fact that L/A + L/n & 1. Until the 
mechanism that determines the location of the sound 
source can be determined, it is not obvious whether 
these results are proof of the feedback tuning mechan- 
ism or merely a result of some other mechanism that 
relates all these characteristic dimensions to the shock 
cell size. 

Analyses have been made of the natural frequencies 
of vibration of incompressible jets, and these have been 
checked with experiment.‘ The relations give the 
TABLE | 


Stagnation pressure for 
natural vibration 


Stagnation pressure at 
without baffles, psia 


Frequency, cps peak noise, psia 


6,000 65 64 
7.250 59.5 56.5 


8.000 58 51.5 
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wavelength of jet oscillations as a function of the size 
of the jet and the Strouhal number. If the method of 
reference 4 is applied to the present experiments, the 
jet wavelength as a function of stagnation pressure is 
found The calculated 
lengths compare favorably with those measured, sug 


as shown in Fig. 12. wave 
gesting strongly that this mechanism of jet instability 
is important in the jet screech problem. The applica 
tion of this incompressible theory to a case of com- 
pressible flow is clearly incorrect. The Mach number 
must have importance in the case of compressible 
flow, which is not indicated by the incompressible 
theory, yet the velocity itself must be of some im- 
portance since it is needed in forming the Strouhal 
number. Since the experiments were all run at an 
approximately constant speed of sound, the distinction 
between Mach number and velocity cannot be de- 
termined empirically 

Thus far, several empirical and theoretical relations 
have been considered which seem to be important. 


It is interesting to attempt to combine them. The 
relations which seem appropriate are as follows: 
n = n(St, d) 
or, for a compressible jet, 
n = n(St, M, d) (1) 
(//n) + (//A) = Co + int. (2) 
An =a (3) 
9= CL, dX of ee (C3 + int.)L (4) 
L= 14, @ (5) 


Since by (3) a a/X, and since, at a given /, a is 


characteristic of all the velocities, (1) becomes 
n = m(A, M, d) 
and hence relation (2) can be written 


l = 1A, M, d, Cy + int.) 





Fic. 10. Shadowgraph of two-dimensional nozzle, Po = 54 
psia: (a) without sound-absorbing material on flange; (b) with 
sound-absorbing material on flange. 
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Combining relations (4) and (5) leads to similar results: 


n = m(A, M, d) 
l = 15(X, M, d, C3 + int.) 


At the stagnation pressures which give maximum 
oscillations, the values of 7 and / calculated as m, /; and 
nm, /, must agree. Three maxima were observed in 
the range studied, so that agreement must occur at 
the three respective pressures. The principal wave 
systems all occur at the same value of the integers in 
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Fic. 11. Output of microphone placed above baffles, with 
speaker baffle arrangement in Fig. 8, as a function of P») and 
speaker frequency. 


relations (2) and (4), so that changes in these integers 
do not play a part in the tuning mechanism. One 
explanation of the jet behavior would be to assume 
that the tuning mechanism depends on the corre- 
spondence of these two functional relations for y and /. 
It should be noted that since the locations of the peak 
amplitudes change with the conditions at the base of 
the jet, these relations must be sensitive to the degree 
of excitation. Another point of view might suggest 
that since the functional relations agree at three points 
they are probably identical, and hence the tuning 
mechanism is not provided by the agreement of these 
relations. 

The five relations presented are based on different 
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degrees of understanding, and not all are equally 
firmly established. The firmest are (3) and (5); 
(3) is almost trivial while (5) is firmly established at 
least for the inviscid core; (1) and (2) have a quasi- 
theoretical basis and have been checked empirically. 
Relation (4) is purely empirical and the basic mecha- 
nism involved is not understood; it might be merely a 
consequence of the other relations, but this seems 
very unlikely since the role of the shock waves would 
then be eliminated. Relation (4) must be con- 
nected with shock-wave effects, although the mecha- 
nism involved here is not clear. 


Conclusions 


(1) The oscillations of an overpressure plane two- 
dimensional sonic jet and the associated sound field 
are directly related to the size of the shock-wave cells. 


?- 

5 st 

: CALCULATEC 

2 - “ MEASURED 

d 

4 et 

“ 
B+ 
230 4 50 60 70 
Po psia 

Fic. 12. Comparison of measured and predicted jet wavelengths 


at various Pp 


(2) Peak amplitudes of oscillation and sound occur 
at particular values of stagnation pressure. 

(3) The interaction of the sound waves with the 
jet near the base are an important part of the instability. 
Shielding the base of the jet from the sound waves 
stabilizes the flow. 

(4) The oscillations occur at frequencies of the 
order of the expected natural jet frequencies. 
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Velocity and Range Considerations for the 
Attainment of Given Intercept Angles 


WAYNE H. TEMPELMAN* 


Convatr- Astronautics, A Division of General Dynamics Corp. 


Summary 


Based on simplified two-body assumptions, a general analysis 
is presented for ballistic trajectories which result in a specified 
intercept angle at a given altitude. The minimum velocity 
required is derived starting with zero velocity and from a circu- 
lar orbit. The range is obtained for initial intercepts at the 
specified altitude. The 
a 150-nm circular orbit with final altitudes of 450,000 ft and 


velocity requirements starting from 


1,028 nm are presented in graphical form 


Symbols 
a = constant defined in Eq. (4) 
b = constant defined in Eq. (4) 
é = eccentricity 
} = radial distance from center of attracting body 
R = dimensionless distance ratio, 7/77 
v = velocity 
J = dimensionless velocity ratio, v/v, 
AV = velocity required starting from a circular orbit 
a = dimensionless radial velocity, V» sin + 
B = dimensionless horizontal velocity, Vp cos 4 
+ = flight path angle with respect to the horizontal 
6 = velocity (AV) angle with respect to the horizontal 
60 = range central angle 
d = Lagrange multiplier 
@ = true anomaly 
Superscripts 
¢ = circular 
d = downward 
i = inward 
I = intercept 
o = outward 
p = perigee 
“u = upward 
Q = initial conditions 


Introduction 


_ ARTICLE IS CONCERNED PRIMARILY with that 
class of trajectories which minimize the velocity 
required for a given intercept angle and altitude or, 
conversely, which maximize the intercept angle for a 
given velocity and altitude. 

A particular application would be the return, with 
small guidance errors, of a data capsule from a satellite 
orbit. Analytical methods are developed in references 
1 and 2 which define the maximum re-entry angle in 
terms of the heating and deceleration limits of the 
capsule during atmospheric descent. As the re-entry 
dispersion is generally minimized for maximized re- 

Received May 3, 1960. Revised and received July 1, 1960. 
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entry angle, the formulas herein presented would be 
useful in selecting optimum orbital ejection conditions 
for specified re-entry angles. 

The usual assumption is made that the system con- 
sists of two bodies, one the attracting body considered 
as a point mass and the other the orbiting body with 
negligible mass as compared with that of the attracting 
body. The only force assumed to act is gravity, hence 
the influence of other bodies, the possible oblateness 
of the attracting body, and atmospheric effects are 


neglected. 


Velocity Requirements 


The velocity of a body in an orbit can be specified 
by a magnitude and a flight-path angle, which is taken 
herein as the angle with respect to the local horizon. 
If, at the initial point, the dimensionless velocity is rep- 
resented by the components a@ (radial) and @ (hori- 
zontal), the initial velocity required starting with zero 


velocity can be expressed as in Eqs. (1). The required 


velocity starting from a circular orbit (a = 0, 8 = 1) is 
specified by Eqs. (2) (see Fig. 1). 
Vo = (a? + B?)'/? Yo = tan! (a/B) (1) 


AV = [a? + (1 — B)?]'? 6 = tan [a/(1 — B)] (2) 


The flight-path angle at the desired altitude will be 
referred to as the intercept angle (y,), and is derived 
from simplified two-body assumptions: 


y, = tan~'[(—R°e? + 2R — 2+ a® + B*)'/?/(RB)] 
(3) 
where R = ”%/r; = (initial distance) /(intercept dis- 


tance). Solving for the radial component of velocity, 


a = (aB? + b)'/? (4) 
where 
a = (R?*/cos? y;) — 1 
b = 2(1 — R) 


Through use of Eqs. (1), (2), and (4), 6 can be ex- 
pressed as a function of either yo or 6: 


B = [b/(tan? y. — a)]'” (5) 


— (6 
tan?6 — a } 








682 JOURNAL OF THE AEROSPACE 
Ve ~ 
Nn ' a OUTWARD INTERCEPT 
/ . ALTITUDE 
“J, \ 


—~_ \ 
> INITIAL ALTITUDE 

| - r ‘, 

| ~ e ° ‘ Vo 


% pig INWARD INTERCEPT 
a as 
\ Ae ALTITUDE 
\ | va 
\ “ 
\ 
\ F 
- wa 
6, 
6 2 , &, 
Fic. 1. Illustration of perigee and range angles. 


The plus sign applies for 6 > 90°, the minus sign for 
5 < 90°; and for 6 = 90°, B = 1. 


Minimum Energy Requirements 


On substituting Eq. (4) into Eq. (1) it becomes 
apparent that the velocity V, for a specified R and y, 
is minimum for the minimum permissible 8. For out- 
ward intercepts (R < 1, b> 0), Eq. (5) is minimized 
for yc = 90° (8 = 0), and for inward intercepts (R> 1, 
b < 0) the value is yy = 0° (a = 0). The correspond- 
ing velocities are found by combining Eqs. (1) and (4): 

R<1: Vo = bi? \ - 
R>1: Vy = (—b/a)"?f 1) 

Minimization of the required velocity (AV) starting 

from a circular orbit is accomplished by minimizing 





oy) 


a] 


a eS 


Illustration of velocity nomenclature and intercept 
geometry. 


Fic. 2. 
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Eq. (2) subject to the constraint on a, Eq. (4). With 
the use of Lagrange multipliers, the following three 
equations are obtained : 

0 = a? — agp’ — b / 

a/{a? + (1 — B)?)"? + 2ha (8) 
0 = (1 — B)/[a@? + (1 — B)?}!? + on,a8) 


I| 


By use of Eqs. (2) and the assumption a # 0), the 
solution is obtained as 


8 = (cos? y;)/R? | 
AVY = 13 — 2R — [(cos? y;) R?]}! 2 (9) 
5 =cos! }[AV? — 2(1 ~ r)V/av}s 


For a = 0, the solution to Eq. (8) is given by Eq. 
(10), with AV and 6 obtained from Eq. (2). 


8 = (-—)b/a)'”? | 
AV =|1 — (—6)/a)'” (10) 
R> 1:6 = 0° R< 1:8 = 180° 

The range of applicability of the solution for a ¥ 0) 
is determined by setting Eq. (4) equal to zero, which 
results in Eqs. (11): 


8 =[1 + (9 — 8R)!/?]/2 
Yta.» = cos "'{R([1 + (9 — SR)"/?]/2)¥f | an 
R> 1: ,< 71 < Yn | 


R< 1: Yq < ¥1< 99° 


The solution minimizing Vy !Eq. (7)] is equivalent 
to that minimizing AV [Eq. (10)| (since AV = 1 — Jj 
when a = y = 6 = QO), both for R > 1.125 and for 
1 < R < 1.125 when jy, is outside the range given in 
Eq. (11) (see Fig. 3). The region where 6 = 0° 
results in a minimum energy intercept originating from 
a circular orbit is shown in Fig. 4. 


Range Considerations 


The above analysis results in those velocity vectors 
which produce intercept at the specified intercept angle 
and altitude. Each velocity will result, in general, 
in intercepts at two different ranges depending on 
whether the velocity is applied upward or downward. 
This is due to the fact that the orbit intercepts a given 
radial distance in two places. 

The intercept point is oriented relative to the initial 
point by the range angle © and to the perigee point by 
the true anomaly ¢ which, in turn, is oriented relative 
to the initial point by the perigee angle @,. There 
are two variables to consider in writing the range equa- 
tions: (1) whether the desired intercept altitude is 
greater or smaller than 7) and (2) the upward or down- 
Fig. 2 shows the 


ward application of velocity is used. 
four possibilities, from which the following ranges of 
angles are inferred for the cases restricted to initial 
intercept, the angles (measured from the last fictitious 
being 


perigee passage in the rotational direction) 
restricted to the principal values: 


Outward Intercepts (r; > ro): 


Upward application of velocity 








Ww 
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\ | 
Box = ? acai 6, 
7 : 80 | —\ i _| 
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Ur a“ . ._ - . 
on velocity components from Eqs. (4) and (1) yields the 
ard. 10 required velocity AV and the range 
ven 
Ok 
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lere 40 i AO 
jua- 
> is The range equations can now be derived: 
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: ieee: : : for the attainment of given intercept angles for a descent from 
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Fic. 6. Velocity requirements starting from a circular orbit for 


the attainment of a given range and intercept angle for a descent 
from 150 nm. to 450,000 ft. 


B = (2/R)' cos y; 

AV = [3 — 2(2/R)"” cos y7}"” (15) 
0 

Omax 27 — cos! [(2 R)(cos? 71) — 1] mas 


cos—! (2 cos? y; — 1) 


For outward intercepts (R < 1) there is a minimum 
6 (maintaining the original orbital direction) that cor- 
responds to a vertical trajectory (yo = 90°, 8B = 0) 
that just reaches the desired altitude. The required 
velocity Vo is given by Eq. (7), and yields 6,in and AV 
from Eqs. (2): 


bmin = tan—'[2(R — 1)]*/?t (16) 
AV = (2R — 1)!2 f aie 


The maximum applicable angle 6 and corresponding 
ranges for upward and downward application of ve- 
locity are determined from the geometry of a straight- 
line transfer: 


Smar = COS 1[(—cos 1) ‘R| l 


: 17 
Ou,a = ¥1 + cos'[(cos y;)/R] J (7) 


where y; > cos! R. 


Figs. 5 to 8 show a typical example of both an out- 
ward and inward intercept. The dashed lines in Figs. 





1,000 av (FT/SEC) 10,000 100,000 

Fic. 7. Possible velocity vectors starting from a circular orbit 

for the attainment of given intercept angles for an ascent from 
150 nm. to 1,028 nm. 
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6 and 8 separate the upward from the downward appli- 
cations of velocity for each y; curve, with the dashed 
lines in Figs. 5 and 7 representing the locus of minimum 
energy intercepts. 


Conclusions 


(1) For inward intercepts with R > 1.125, the mini- 
mum energy-transfer trajectory starting with zero 
velocity or from a circular orbit is obtained by a hori- 
zontal application of velocity. 

(2) For inward intercepts with R < 1.125, there are 
two ranges of intercept angles, one starting from zero 
degrees and the other ending at 90 degrees, where the 
first conclusion holds true. When the intercept angle 
lies between these two ranges, minimum energy trans- 
fers starting from a point with zero velocity are ob- 
tained with a horizontal application of velocity, whereas 


320 SS 
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Fic. 8. Velocity requirements starting from a circular orbit for 


the attainment of a given range and intercept angle for an ascent 
from 150 nm. to 1,028 nm. 


in starting from a circular orbit minimum energy 
transfers are obtained with a nonhorizontal applica- 
tion of velocity. 

(3) For outward intercepts, minimum energy trans- 
fers starting with zero velocity are obtained by a verti- 
cal application of velocity, resulting in a range of zero 
degrees. Minimum energy transfers starting from a 
circular orbit are obtained by a horizontal forward 
application of velocity (6 = 180°) for a limited range 
of intercept angles starting at zero degrees. Outside 
of this range the angle of velocity application decreases 
to some angle less than 90 degrees. 

(4) Inward intercepts for a given y, result in a limited 
range coverage, range angles near zero degrees and 
near 360 degrees, being excluded. 

(5) Outward intercepts, for a range of intercept 
angles starting at zero degrees, allow complete coverage 


(Continued on page 701) 
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Aerodynamic Noise in Supersonic Wind 
Tunnels' 


JOHN LAUFER* 
California Institute of Technology 


Summary 


Hot-wire measurements in the free stream of a supersonic 
wind tunnel were made in the Mach number range of 1.6 to 5.0. 
It is shown that the mass-flow fluctuations increase very rapidly 
If the fluctuation field is as- 

an assumption that is con- 


with increasing Mach number. 

sumed to consist of sound waves 
sistent with the measurements—the sound intensity is approxi- 
mately proportional to M+, within the range of the experiments. 
Furthermore, the orientation of the field is found to be different 
from the Mach line direction; it corresponds to a sound-source 
velocity of approximately one-half the free-stream velocity for 
the higher Mach numbers. It that the turbulent 
boundary layer along the nozzle and the tunnel walls is respon- 


is shown 


sible for this sound field. 


Symbols 
a = sound velocity 
e = voltage across hot wire 
Aer = sensitivity coefficient for total temperature 
Aem = sensitivity coefficient for mass flow 
m = mass flow 
M = Mach number 
p = pressure 
r = Ae,,/Aep 
R = correlation coefficient 
Re/in. = Reynolds number per in. based on free-stream 
conditions 
T = static temperature 
Tr = total temperature 
u = velocity in the mean flow direction 
u,, = particle velocity normal to the wave front 
i, = source velocity 


Y- 1 1 
a = 1+ ee ur) 


8 = a(y — 1)M? 

Y = ratio of specific heats 

6* = boundary-layer displacement thickness 

0 = angle between the normal to the wave front and the 
direction of flow 

p = density 

( = instantaneous values 

() = mean values 

¢) = root-mean-square values 
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Introduction 


is IS WELL KNOWN that in a low-speed wind tunnel 
certain types of experiments are strongly affected 
by the turbulence level of the free stream. The classi- 
cal example is the drag measurement of a sphere. 
Prandtl-Wieselsberger' pointed out that the differ- 
ence in drag values obtained in various wind tunnels is 
due to turbulence in the free stream. Later, Taylor® 
obtained a parameter correlating the measurements. 
Dryden and co-workers’ carried out systematic experi- 
ments and substantiated Taylor’s theoretical work. 
When a better understanding of turbulence produced 
behind grids had been gained, effective methods were 
devised to reduce the turbulence level in the free 
stream.* 

In all the preceding cases the disturbances were pre- 
dominantly velocity fluctuations having vorticity. 
Since no cooling systems were required, very few tem- 
perature fluctuations were present in the stream. 
Furthermore, disturbances due to noise could be mini- 
mized by proper design of the propeller and its shaft. 

It soon became apparent that in supersonic wind 
tunnels the fluctuations differed in source and nature 
from those experienced in low-speed tunnels. In 1954, 
the author noted® that large velocity fluctuations in 
the supply section do not influence the transition Reyn- 
olds number in the free stream at high Mach num- 
bers. He suggested, as a possibility, that sound-type 
disturbances originating from the turbulent shear layer 
on the tunnel walls might be present in the free stream. 
Later, in 1956, some wind-tunnel measurements at the 
NACA Lewis Laboratory indicated that very large 
fluctuations are present in an M = 3 flow.t The most 
comprehensive treatment of the subject is given by 
Morkovin,® 7 who discusses all the possible sources of 
disturbances in a supersonic wind tunnel, points out 
the difficulties of measuring these fluctuations, and 
gives an approximate value for the turbulence level 
encountered in the wind tunnel used in his investiga- 
tions. The present investigation was undertaken in 
connection with experiments relating to the stability 
problem of the supersonic laminar boundary layer. 
From experience obtained in low-speed work, it was 
known that the detection of the unstable Tollmien- 
Schlichting oscillations is difficult, if not impossible, 
when the free-stream disturbance level is high. For 





t Private communication from J. C. Evvard, Lewis Labora- 
tory, Cleveland, Ohio. 
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Fic. 1. Schematic diagram of the tunnel. 
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Fic. 2. Variation with Mach number of boundary-layer tliick- 
ness on the tunnel wall. 
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Fic. 3. Boundary-layer thickness distribution in a tunnel cross 
section. 


this reason, considerable time and effort were spent 
to study the fluctuations in the free stream and to find 
means of reducing their amplitude. 


Equipment 


Wind Tunnel 


Most of the experiments were carried out in the Jet 
Propulsion Laboratory 18 XX 20-in. supersonic wind 
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tunnel. The peak Mach number of the tunnel is 5. 
Dimensions and details of the supply section, nozzle, 
and test section are shown in Fig. 1. Screens | and 2 
in the cone section are used to prevent separation. 
The air filter is a commercially available Airmat paper 
filter} with a drag of about 50 to 100 times the dynamic 
pressure. The first four damping screens are 20 mesh 
with ().016-in. diameter steel wires; the last two screens 
are 30 mesh with 0.003-in. diameter wires. In the 
supply section all corners are faired. 

The flexible-nozzle walls are of the continuous third- 
derivative design.* The mean pressure variation in the 
working section is less than 2% for the entire Mach 
number range. Details of the mean flow distribution 
at various Mach numbers are given in reference 9. 
The boundary layers on the tunnel walls at the throat 
and downstream from the throat are turbulent under 
all flow conditions except when otherwise stated. Fig. 
2 gives boundary-layer displacement thicknesses meas- 
ured on the center of the upper flexible wall as a func- 
tion of Mach number. The approximate shape of the 
thickness distribution around the circumference of the 
tunnel test section is shown in Fig. 3. 

The turbulence level in the supply section was found 
to be due to velocity fluctuations only, no temperature 
fluctuation being detected. The turbulence level was 
1 per cent for all Reynolds numbers and Mach numbers, 
except at M/ 2 4.5, where fi/a# = 0.5 percent. 


Hot-Wire Probe 

All the hot wires employed in the experiments were 
made of 90 per cent platinum and 10 per cent rhodium 
with diameters of 0.00005 or 0.0001 in. and correspond- 
ing nominal lengths of 0.015 or 0.020 in. They were 
mounted on a probe as shown in Fig. 4. 
body of the probe was hollow and diamond-shaped 
with a 5-degree wedge attached. A thin razor blade 
0.060 in. wide and 0.004 in. thick was soldered to this 
wedge. The razor blade was coated with an electrically 
insulating varnish; two jeweler’s broaches were glued 
to the blade, and a final coat of varnish was baked over 
the surface. The broaches extended 0.050 in. in front 
of the thin, razor-sharp leading edge and were spaced 
either 0.015 or 0.020 in. apart. Lead wires from the 
broaches were fastened to the rear of the wedge and 
passed through the center of the diamond-shaped piece 
to the gold-plated pins which plugged into a vertical 
With this mechanism 


The main 


traverse mechanism (Fig. 5). 
the hot wire could be placed close to the upper wall of 
the tunnel in order to protect the wire from the rough 
flow conditions which obtain during the starting of the 
wind tunnel. With this traverse, the wire could also 
be placed at any vertical position in the tunnel. 


Hot-Wire Amplifier 

It is more difficult to meet the design requirements 
of the hot-wire equipment used for fluctuation measure- 
ments in supersonic flows than in low-speed flows. 


+ Manufactured by American Air-Filter Co., Louisville, Ky 
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Much higher frequencies are expected in supersonic 
flows, and an amplifier with a frequency response of 
several hundred kilocycles is therefore needed. Low- 
noise input circuitry is especially critical in such an 
amplifier because of the greater frequency bandwidth 
and lower level of the fluctuations expected. The com- 
mercially available hot-wire set* satisfies these re- 
quirements and has been used throughout the present 
experiment. 


Experimental Results 


Early Experiments 


It is well known that the transition Reynolds num- 
ber on a model placed in a wind tunnel is sensitive to 
the turbulence level of the tunnel free stream. This 
fact has been utilized by the author® to obtain some 
qualitative concept of how the free-stream turbulence 
level changes as the tunnel Mach number increases. 
Of course, any conclusion arrived at from this rela- 
tively crude method has to be weighed very carefully, 
since the stability of the boundary layer is also a 
function of the Mach number. However, although 
the available theoretical calculations” indicated that 
in the low supersonic range the boundary layer be- 
comes more stable as the free-stream Mach number 
increases, the experiments indicated a sharp decrease 
in the value of the transition Reynolds number. Thus, 
it was tentatively concluded that the turbulence level 
of a supersonic tunnel increases with increasing Mach 
number. And it was further concluded that the tur- 
bulence is not due to fluctuations originating in the 
supply section, since a tenfold increase in the velocity 
fluctuation there did not affect the transition Reynolds 
number of the model at WM > 2.5. The velocity-fluc- 
tuation increases did, however, decrease the transition 
Reynolds number appreciably at lower Mach numbers. 


Hot-Wire Measurements 

To gain a better understanding of the free-stream 
fluctuations, a more direct method of measurement is 
required. Thus far, the most successful is the hot-wire 
technique, with which the magnitudes of certain types 
of fluctuations can be obtained. The main difficulty 
with this technique is that these fluctuations represent 
an algebraic combination of fluctuations of several 
basic physical quantities, and, with the exception of 
some special cases, it is not possible to calculate the 
magnitude of the fluctuation of each basic quantity 
separately. Specifically. the hot wire in a supersonic 
stream responds to both mass-flow and total-tempera- 
ture fluctuations, and the instantaneous voltage across 
the wire may be expressed!! 


e’ = Aer(Tr'/Tr) — Aem(m'/m) (1) 
where 


Tr'/Tr = af T'/T) + B(u'/a) 


* Model 50B, manufactured by Shapiro and Edwards, South 
Pasadena, Calif 





Fic. 4. Hot-wire probe 
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Fic. 5. Hot wire and total head probe 
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Fic. 6. Mode diagram for WV = 1.6 


and m'/m (u’/i) + (p’/p) 


and Ae, and Aer are sensitivity coefficients that are 
determined by mean flow conditions and the mean 
temperature of the hot wire. The measured mean- 
square value of the voltage may be written in the form 
e”? T;"? mT; ,m” 


-= — — 2rRar —p + 1 (2) 
1 L 


ni; m? 
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/ 
where r = Ae,/ Aer 
I Jas 
and Rar = m'T'7/mT 7 


If under a given flow condition e’? is measured for three 
different values of r (i.e., at three different mean wire 
temperatures), the values of 7, m, and Rnrp may be 
calculated. In practice, the measurements are carried 
out at more than three wire temperatures, and from 
plots of @/Aer vs. r the values of the three unknowns 
may be obtained more accurately by a fairing tech- 
Figs. 6 through 10 show some typical results 
The calcu- 


nique. 
obtained in the 18 X 20-in. wind tunnel. 
lated mass-flow fluctuations are plotted in Fig. 11. 
Their large magnitudes at the high Mach numbers are 
quite startling: they are approximately 50 times 
higher than those measured in a low-turbulence, low- 
speed wind tunnel. 

One remarkable feature of these fluctuations is that 
without exception the correlation coefficient R,,7 is 
found to be approximately —1 (the mode diagrams in 
Figs. 6-10 are always straight lines within the accuracy 
of the measurements)—that is to say, the mass-flow 
and total-temperature fluctuations exhibit a high de- 
gree of correlation. This is a strong restrictive condi- 
Thus, the fluctuating field in the free stream must 
Unfortunately, from the 


tion. 
be of a very special kind. 
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hot-wire data alone it is not possible to obtain any 


more definitive information. What can be done, 
however, is to postulate a certain simple type of fluctu- 
ation field and to determine whether the calculated 
hot-wire response (mode diagram) is consistent with 
the measured response. It is clear that the field can- 
not be due to pure velocity fluctuations associated 


with vorticity. In this case, 


m'/m = u'/tt 
and Tr'/Tr = Bu'/i 
Consequently, 
Rap = m'Te'/mMTr = u"?/u"? = +1 

This expression, in contrast with the results, corre- 
sponds to mode diagrams having straight lines with 
negative slopes. 

Assuming pure temperature fluctuations that are 
not associated with pressure changes, i.e., entropy fluc- 
tuations, then 
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This gives the proper sign for the correlation coefficient 
R,,,- However, rewriting the hot-wire equation in 
terms of the r.m.s. temperature fluctuations 


@/Aer = (a + r)T/T 


gives r = —a for é/Aer = 0. This condition is not 
satisfied by the measurements. 

As a matter of fact, the existence of these two types 
of fluctuations or their combinations can be ruled out 
within the framework of the linear theory'” using simple 
physical arguments. Neither in the nozzle nor in the 
working section does the tunnel free stream contain 
turbulent shear fields or shock waves. Therefore. 


since both -vorticity and entropy fluctuations are 
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convected along streamlines, they could only exist in 
the working section through convection from the supply 
section. However, as noted earlier, the temperature 
(or entropy) fluctuations are undetectable here; thus, 
they must be negligibly small in the working sec- 
tion also. Since the per cent of velocity fluctuations 
associated with vorticity in the supply section is greatly 
diminished (to less than 0.01 per cent) through the 
nozzle because of the large ratio (40 at JJ = 1.6 and 
1,500 at M = 5.0) of mean velocities in the working 
and supply sections (see the calculations of Tucker in 
reference 13), they cannot be the source of the high 
fluctuation level measured in the working section. 
The only other simple fluctuating field'* is a pure 
sound field in which the isentropic relation between 


the fluctuating quantities holds. Thus, 
p' p’ . 2 
—=y- = - = 
P - yore 


and 

mou’ 1 p’ T;’ , (2 “) 
ee: =ey7-1) (2+ 

= i a T, a(y ) a = 2 


The hot-wire equation now becomes 
, / , , 
u b u 
e’ = a(y — 1) (2 + M? -— ) Aer — (é — -- — i Ae. 
YP u yD ou 


, 


é 
or — = (8 — r) 


u 
Aer u 


+ [a(y -— 1) - 7] a 
YP 


The mean-square voltage is 


9 9 


e’2 
Aer” 
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number. 
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It is not difficult to show algebraically that the condition R,,7 = —1 corresponds to R,, = —1; accordingly, the 


final result is 


il p 
= (8 — r) te laly ~ 1) — #] ~ 


é 
Aer 


The assumption of a pure sound field with waves 
moving against the mean stream (X,, = —1) can thus 
be made consistent with the measurements. Using 
the mode diagrams obtained from the measurements, 
it is now possible to calculate separately all the fluctu- 
ating quantities present in the free stream. Fig. 12 
shows the r.m.s. pressure fluctuations (made dimen- 
sionless by the dynamic pressure) as a function of the 
free-stream Mach number. It is to be noted that 
the mean pressure level of the tunnel (that is, the value 
of the Re/in.) influences the free-stream fluctuation 
amplitudes. 

Providing that the preceding interpretation of the 
measurements is correct, further information can be 
extracted from the hot-wire data. It can be shown 
that the sound has a preferred orientation which is not 
in the direction of the characteristics; in other words, 
the measured pressure variations or at least their major 


hi p p ii 
—1)(Ww-- _ 
cv o(wr-F)+(B-2)r 





portion are not caused by fluctuating Mach waves. 
Since the average wavelength of the field can be esti- 
mated to be of the order of an inch or less, it is safe 
to assume that the measurements are made in the far 
field of the sound source, and the waves may be as- 
sumed to be plane. 

For a plane wave, the relation between u,,’, the par- 
ticle velocity normal to the wave front, and p’ may be 
written 

| wn’| 1 p’ 
i Mvyp 
However, it should be remembered that the hot wire 
is sensitive only to velocity variations in the direction 
of the mean flow such that 


| | 
u’ = |uU,n’| cos 8 














690 JOURNAL OF THE AEROSPACE 
10 T T ] Ci a 
| | | 
| | } 
al. @ Fe/in = 340,000| | | i 2 
o Re/in = 90,000 
© | | ° ° | 
& eam tele Fo on, 
BS | ? | | | 
P 4r + & — = : ¢ +— a 
Iq | 
| 
2} | { 4 H = 
) ae | | _ 
2 3 4 5 
M 
Fic. 12. Variation of pressure fluctuations with tunnel Mach 


number. 


where @ is the angle between the normal to the wave 
front and the direction of the mean fiow. Therefore, 


u’/n = (1/M)(p'/yp) cos 6 (5) 
For a Mach wave, cos 6 = —1/.M], and the r.m.s. pres- 
sure fiuctuation due to a fluctuating Mach wave is 
given (see also reference | 1) 


b/ yp = M*(ii/it) (€) 
Examining Eq. (4), 
e/Aer = |(M? — 1)/M*\(p/yp)r 


Thus, in the mode diagram the straight line goes 
through the origin. This, however, is not the case for 
the present measurements. 

In order to determine whether the accuracy of the 
measurements is sufficiently high to rule out the exis- 
tence of fluctuating Mach waves in the free stream, the 
following experiment was performed. A tape 0.004 in. 
thick was placed across one of the nozzle walls of the 
wind tunnel at such a position that the weak shock 
wave originating from the tape crossed the hot wire. 
The shocklet was, of course, fluctuating because of the 
presence of a turbulent boundary layer on the nozzle 
wall. If the velocity and pressure due to these fluctu- 
ations are designated by wm,’ and p,’, then, by Eq. (6), 


rx 10° 


&/ae 











Fic. 13. Mode diagram for a fluctuating Mach wave superim- 
posed on free-stream fluctuations. 
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p»" yp = M?(1u2’ nu) 


Furthermore, it is easily shown that in the resulting 
mode diagram the intercept at r = () should have the 
saine value as measured in the clean tunnel without the 
fluctuating shocklet. That is, writing e. for the sum 


by = ea, 


| = tana (ar oe a ) - ( i. ) 
Aer _|, =0 TT p Aer], =0 


where the subscript 1 corresponds to the clean tunnel 
condition. 

The mode diagram in Fig. 13 shows that at r = (0) 
the value of @/Aer did not change with the shocklet 
present. 

This experiment supports the reliability of the hot- 
wire technique and indicates that sound fluctuations 
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the orientation of which is not in the direction of the 
characteristics are indeed present, which means that 
the source of this field cannot be fixed but must be 
moving with respect to the tunnel reference. 

Using Eq. (5), and remembering that R,, = —1, 
it is possible now to calculate the most probable orien- 
tation of the waves and, in turn, to obtain a mean 


propagation velocity of the sound source u,;. Since 
(i — us)/a = —1/(cos 6) 
ih. l b/yp 
“Teele ai ih 
ut M cos 6 M*ii/ tt 


Fig. 14 represents this ratio measured at various free- 
stream Mach numbers. No values are shown at lower 
Mach numbers; here the determination of #,/a# be- 
comes less accurate because of the relatively lower 
amplitudes of the fluctuations. It should be mentioned 
here that the computed value of “, is rather sensitive 
to the assumption that a// the sound fluctuations are 
emitted from the same source. If, for instance, some 


Mach-wave fluctuations are also present—a condition 
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that cannot be completely ruled out, especially in the 
low Mach number flows—the computed value of 
a, would increase.* The points in Fig. 14 thus repre- 
sent a lower limit for the propagation velocities. 


Discussion 


The preceding interpretation of the hot-wire meas- 
urements shows the fluctuating field existing in the 
wind tunnel to be a sound field the intensity of which 
increases rapidly as the tunnel Mach number is in- 
creased. (Incidentally, the intensity was found to be 
uniform across the tunnel—less than | per cent vari- 
ation—especially for the higher Mach number flows: 
as the flow Mach number decreased, the extent of the 
uniform field decreased toward the center position of 
the tunnel.) Furthermore, the field is not isotropic, 
but has a preferred direction. From this direction 
the source velocity was computed and was found to be 
between ().4 and 0.5 times the free-stream velocity. 

The question then arises as to what mechanism can 
produce such a fluctuation field. It has already been 
shown that the field does not originate from the supply 
section, where the sound fluctuations are negligibly 
small. Mechanical vibrations are also ruled out, since 
the frequencies encountered are very high. It is pro- 
posed, therefore, that the reason for this intense sound 
field is the presence of turbulent boundary layers along 
the tunnel walls: the fluctuations are identified with 
the aerodynamic noise produced by the turbulent 
boundary layers on the walls. In the following dis- 
cussions arguments and evidences will be described 
which support this contention. 

(a) It may be argued that if the boundary layers on 
the tunnel walls are indeed responsible for the fluctua- 
tions picked up by the hot wire, and if the hot wire is 
“protected’’ from the tunnel walls by an apprc- 
priately placed shield, a reduced voltage fluctuation 
should appear across the wire. The following experi- 
ment was therefore performed. A flat plate was placed 
on the center line of the tunnel spanning the working 
section from wall to wall. A hot wire was moved in 
the free stream over the plate in such a way that it 
could traverse ahead of the plate leading edge and 
behind the leading-edge shock wave. In the latter 
position the wire was shielded by the plate from the 
radiation coming from the boundary layer on the bot- 
tom wall. Of course, it was still exposed to the re- 
flected field of that layer and to the direct field of the 
layers on three other walls. Nevertheless, Fig. 15 
clearly shows the mean-square voltage fluctuation 
across the wire (which is roughly proportional to p*p*) 
decreased by about 20 per cent. (The wire sensitivity 
change was negligible ahead of and behind the shock.) 
The large peak in Fig. 15 corresponds to the position 
at which the wire crossed the leading-edge shock. 
Ahead of that position, e* corresponds to the free- 


* The author is grateful to Dr. M. V. Morkovin who pointed 


out this to him 
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Fic 15. Free-stream fluctuations with a ‘‘shield.”’ 


stream fluctuation level; behind it, the reduced values 
of @* indicate the shielding effect of the plate. 

(b) Adopting the suggestion of Liepmann' that the 
sound field due to a turbulent boundary layer may be 
considered the result of the displacement thickness 
fluctuation of the layer, it is expected that the sound 
intensity increases with increasing mean boundary- 
layer thickness. Indeed, the experiments show that 
for a constant Mach number fiow with an increase in 
Reynolds number per inch of the tunnel, which corre- 
sponds to thinning of the wall boundary layers, the 
sound-field intensity decreases (Fig. 12). 

(c) The recent experimental work of Willmarth” 
has indicated that the pressure fluctuations on the 
wall due to a turbulent boundary layer have a ckarac- 
teristic propagation velocity which is 0.8 and is con- 
stant with Mach number in the range of the experi- 
ments (up to 1J = 0.67). Although the connection 
between the wall and free-stream pressure fluctuations 
produced by a turbulent layer is not understood as 
yet, it is clear that they must be interrelated. Con- 
sidering Willmarth’s results, it is not surprising, there- 
fore, that in the present work it is possible to identify 
a source velocity which is different from the free- 
stream Mach number. 





Fic 16. Shadowgraph of model in free flight 
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(d) Schlieren photographs obtained in a wind tunnel 
cannot reveal the pressure radiation field in the free 
stream because of the presence of turbulent boundary 
layers on the tunnel windows. However, some 
schlieren photographs taken in free flight'® do reveal 
that it is possible to recognize an average orientation 
of the sound field emanating from the turbulent 
boundary layer which does not coincide with that of 
the Mach lines (Fig. 16).* As a matter of fact, from 
a rough estimate of sound-field orientation, a source 
speed is obtained which is about half the free-stream 
velocity and is in agreement with the present findings. 

(e) Finally, it was possible to operate the wind tunnel 
at pressures low enough (Re/in. = 26,000; M7 = 4.5) 
that the boundary layers on the four walls in the test 
section were observed to be laminar. Under these con- 
ditions, the mass-flow fluctuations decreased by a 
factor of 10 (to about 0.11 per cent). The mode 
diagram indicated furthermore that these fluctuations 
were primarily temperature disturbances in the flow. 
This experiment thus proved definitely that the source 
of the large free-stream turbulence is indeed the turbu- 
lent wall boundary layers. 


Conclusions 


Hot-wire measurements made in the working section 
of a supersonic wind tunnel showed that the ampli- 
tudes of the mass-flow fluctuations increase very 
rapidly as the flow Mach number is increased. They 
further indicated that the fluctuation field is of a very 
special kind: one in which the mass-flow and _total- 
temperature fluctuations show a high degree of cor- 
relation. 

These measurements are interpreted to be the result 
of a sound field which is proportional in intensity to 
MM‘, within the range of the experiments, and in which 
the sound waves have a preferred orientation. A 
source velocity can be calculated which was found to 
be approximately one half the free-stream velocity at 
the higher Mach numbers. 

Finally, it is proposed and proved that the sound 
field is produced by the turbulent boundary layers 
on the tunnel walls. 

The nature of the free-stream turbulence described 
here is believed to be characteristic of all supersonic 
wind tunnels that have turbulent boundary layers 
on the walls and good flow conditions in the nozzle and 


* Fig. 16 originally appeared as Fig. 4c of NACA Technical 
Note 4235 (Observations of Turbulent-Burst Geometry and Growth 
in Supersonic Flow by C. S. James). The author is grateful to 
the National Aeronautics and Space Administration for per- 
mission to republish it here. It is to be noted in this figure that 
the wavelets near and along the leading-edge shock are due to 
turbulent spots and are to be disregarded; instead, reference 
should be made to the fine wave structure in the vicinity of the 


turbulent boundary layers. 
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supply section. Some measurements were carried 


out, for example, in a smaller (12 X 12-in.) tunnel 
with similar results. 

The presence of the sound field in the wind tunnel 
may be utilized to great advantage in the study of the 
aerodynamic noise produced by supersonic turbulent 
boundary layers. Thus far such work has not been 
done experimentally. 

On the other hand, the presence of the sound field 
handicaps any experimental work on stability and 
transition problems which arise at high Mach numbers. 
On the basis of the present results, any measurements 
of transition Reynolds number obtained in wind tunnels 
with turbulent wall boundary layers should be accepted 
with great reserve. 
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Transter From an Arbitrary Initial Flight 
Condition to a Point Target 


HARVEY L. ROTH* 
Bell Aerosystems Company 


Summary 


In Part I of this paper is developed a noniterative method for 
obtaining the minimum impulse necessary to reach a stationary 
point target from an initial flight condition where the target need 
not be in the original plane of motion. The degenerate cases of 
minimum energy trajectories and departure from a circular orbit 
are given special attention. In Part II are obtained several sets 
of partial derivatives essential to an evaluation of range and/or 
flight-time errors at target. Special consideration is given to the 
singular points at which the Jacobian determinant associated 
with the transformation from one set of derivatives to another 
vanishes. Although the two parts constitute separate problems 
both are needed for an overall evaluation of particular missions, 


PART I 


Minimization of Fuel Required to Reach a Point 
Target From Any Initial Flight Condition 


(A) Preliminary Remarks 


It is important to find the most economical manner 
in which to move a vehicle from one point to another. 
This paper considers the single-impulse transfer from 
an arbitrary initial flight condition to a specified point 
In the first part, the minimum impulse 
transfer is obtained. In Part II the pertinent error 
derivatives and their singularities are discussed. The 
entire analysis deals with motion in an inverse-square 


destination. 


force field. 


(B) The Minimum-Energy Trajectory 

It is well known that a spherical body sets up a gravi- 
tational field where the force F (per unit mass) is of the 
following form (reference 1, pp. 95-97): 


F = F(r) = K/r’ (1) 


where A = constant for a given sphere and ¢ is the 
geocentric radius or distance from the center of the 
sphere. It is further known that a vehicle under the 
influence of only the above force will move along a 
conic section, i.e., a curve described (reference 1, pp. 
227-229; reference 2, pp. 1A1-—1A7, etc.) by the fol- 


lowing equation: 


p - p 
1 — ecos 6; 


ri, = (2) 


1 — ecos 0; 
where p is the semi-latus rectum, e is the eccentricity, 
and @ is the angle between the major axis and the geo- 

Received May 3, 1960. Revised and received July 2, 1960. 
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centric radius vector. 
é=0° 
If ’; denotes the vehicle speed at the point 7, it can 


Since p, e > 0, r = Phar When 


be shown that? 


p = 7,8; cos” ¥; | 


: ' (3) 
e= V1 — 82 — B;) cos?y;! ; 


where 8; = V,*r;/X; X = gr? = gore? = const.; g = 
acceleration due to gravity at the surface of the sphere; 
and y = the angle between the velocity vector and the 
local horizontal plane. If the projection of the ve- 
locity V; on the radius vector points away from the 
center of the sphere, y is positive; otherwise 7 is nega- 
tive (see Fig. 1). 

The range RX from point 7 to point f is defined as the 
length of the projection on the sphere’s surface of that 
portion of the trajectory between i and f. If the radius 
of the sphere is denoted by 7, then R = (6; — 6;). 

The range angle @g is defined by the following ex- 
pression : 


Or => R %= 6; = 6; (4) 
and therefore: 


cos 6; = cos (6; + 9p) 
= cos 6; cos 6g — sin 6; sin Oz (5) 


Combining Eqs. (2), (3), and (4) yields the following 


result: 


t The local horizontal plane at a point 7 is the plane perpen- 


dicular to 7; 
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Basic trajectory variables 


Fic. 1. (left). 
velocity vector diagram. 


Fic. 2. (right). Two-dimensional 
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A(1 — cos 6p) 


Vi = ° | 
i n COS Vilri cos yi + ry sin y; sin 0g — r; cos y; Cos On| 
A(1 — cos OR) | 
n cos yi[r; cos y; — 7, cos (6p + y¥;)| 
l — cos 6, 
or Bi = 


cos yi[n cos y; — cos(@z + vi) |] 


where @; is defined in Eq. (3) and n = 7;/ry. 
The minimum speed at point 7 necessary to move a vehicle from point 7 to point f through a specified angle @z (or 
range R) is found by solving either of the following equations for y;: 


dV?7/dy; = V;dV,/dy: = 0, 
or equivalently dp,/dy; = O (7) 


where, in general, |’; ¥ 0. 
From Eqs. (6) and (7) it follows that the minimum velocity is applied at an angle y,, given by the following ex- 


” — cos 0 7 — cos 0p\’ 
tan ¥m = — - ¥( : - + I (8) 
sin Op sin Op 


n 


An examination of derivatives of l’; of various orders with respect to y; shows that V;—> © when d"V; dy;"—> « 
forany”n. Itis further seen that all the derivatives of V’; exist and are finite at all other values of y;. 

However, not all solutions to Eq. (8) need represent physically realizable situations. 

From Eq. (6) it is seen that the points at which l’; > © are given by: 


" COs 8p — 9 
7: = +90 and Yi = arctan : (9) 
sin 6p 


pression : 


The physically impossible situation V’;? < 0 occurs for tan y; < (cos 0g — n)/ (sin Oz). 
It further follows from Eq. (6) that the plot of V’,? versus y; does not go through |’,°> = 0. Therefore if a V;, y; 
combination exists for the 7;, r;, 9g combination in question, then a minimum J’; also exists whose value is determined 


by solving Eq. (5). 


Since the energy / (per unit mass) associated with a trajectory is given by 


a yee VP 
= = = — — = const. 
'r 


bo 
mal 
io 


where r; is specified, minimizing V’; minimizes the energy—i.e., the trajectory chosen above is, in fact, the minimum- 
energy trajectory. 
(C) The In-Plane Transfer From Flight Condition to Point 


Consider a vehicle having a velocity of magnitude J’y) at an angle yp. It is desired to change impulsively the mag- 
nitude and direction of this velocity by adding an increment AV at an angle 6 to the original velocity. The vector AV 
is related to V, and the resulting velocity V, by the following equation (see Fig. 2): 


(AV)? = V2 + Vo? — 2VoVi cos (¥i — Ye) (10) 
sind = (V,/AV) sin (yi — yo) = VB, AB sin (yi — Yo) 
cos 6 = (1/AV)[V; cos (y: — vo) — Vo] = VB; A8 cos (yi — Yo) — V By AB (11) 
where AB = (AV)?r,/xr 
If a vehicle has a velocity Vp at the point 7, then there is a nondenumerable set of values for the AV necessary to 
move the vehicle to a point fin the original plane of motion. Substituting Eq. (6) into Eq. (10) gives Al’ as a fune- 
tion of y; as follows: 


sis A(1 — cos Op) = : A(1 — cos Op) 
(AV)? = + Vo? — 2Vo cos (yi — Yo) 
n cos y;[7; COS Yi — r,cos (0g + 7i)| n COS Y;[7; COS Yi — r,cos (0x + ¥i)} 


or 
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Fic. 3. (left). Spherical range triangle. 
Fic. 4. (right). Three-dimensional velocity vector diagram. 
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Carpet of ‘“‘dimensionless flight time’ ./A/r;* t vs. n 
and 6p for a circular orbit 


052 
0: 
0. 
0.40 
0.36 


0.32 
028 Op DEGREES 










3.0 


Bi 024 
020 

0.16 

012 

008 

004 

0 


Fic. 9. Carpet of 8; vs. » and 6 for a circular orbit. 
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! 
| 


1 — cos Oz | Bo (1 — cos Op) 


Ag => a Bo = 2 cos (¥: -_ Yo) 


(12) 
cos yi[n cos ¥; — cos (On + 7:)] Vccs 79 cos 7, — cos (2 + 01 


As before, the minimum AV is found by setting d(AV’)?/dy; = 0 [see the discussion following Eq. (8)|. Further- 
more, it is shown in reference 3 that minimizing AV results in minimizing the fuel expenditure when Al’ is added 


impulsively. 


(D) The General Three-Dimensional Problem 

In general, an initial velocity Vy will be specified along with the in-plane range Kp and the lateral range L. If 
6p = Ro/r, 6, = L/r, where ro is the radius of the reference sphere, it then follows from spherical trigonometry 
that the total range XK is given by the following expression : 


cos 6p = cos (R/7r) = cos @p cos 6, (13) 


The angle a between the initial plane (the plane of V,)) and the final plane (the plane of 2) is given by (see Fig. 3): 
cot a = sin 6p cot 6, (14) 


It is further noted that the range RX can be achieved by obtaining range FX in the Vp plane and then rotating the 
vector V; such that its magnitude and angle with the local horizontal is unchanged. This is achieved by rotating the 
horizontal projection of V; through an angle a (see Fig. 4). 

In summary, the method is as follows. The proper Vj, y;:,'7;, 7; combination is used in order to obtain a range R 
in the plane of V;. The plane of V; is then rotated through an angle a by adding a velocity V, to V; in such a manner 
that the resultant velocity V; has the same magnitude and y as V;. Since |V;/ sin y; = |V;/ sin y; = |; sin y,, it 
follows that both ends of the vector V, are equidistant from the horizontal plane—1.e., V, is in the horizontal plane. 

From simple trigonometric considerations, the magnitude of V, is given by the following expression: 


V, = 2V; cos y; sin (a/2) (15) 
The total velocity increment AV7 is the vector sum of AV [see Eq. (10) and the preceding remarks] and V,, where, 
from the previous discussion, V, lies in the local horizontal plane and AV lies in the local vertical plane containing the 


vector Vp. 

Consider the following reference system. Let the x-axis be the intersection of the local horizontal and local 
vertical planes containing V;; the y-axis is the perpendicular to the x-axis in the horizontal plane (i.e., the perpendicu- 
lar to the local vertical plane) and the z-axis is a line perpendicular to the above two axes (or the perpendicular to 
the local horizontal plane). 

From the previous discussion, the x, y, components of V, and AV are as follows (see Figs. 2 and 4): 


AV, = AV cos (x — 6 — yw), AV, = 0, AV. = AV sin (x — 6 — wi 


‘ ' ars (16) 
Vicc= V, 666 ; Vy = V, sin , V.. = O 
9 ») 


de ie eS - | 
The magnitude of the total velocity increment AV7 is therefore given by the following expression: 
(AVr)? = (Vez + AV,)? + (AV)? + (Vey)? = [Vz sin (a/2) + AV cos (w — 6 — yo))? + 
[AV sin (wf — 6 — yo)]? + [V. cos (a/2)]? = (AV)? + V,? + 2V, AV sin (a/2) cos (w — 6 — yo) 


(17) 
Substituting Eqs. (10), (11), and (15) into Eq. (17) yields the following result : 
(AV7r)? = Vo? + Vi? — 2VoV; cos (yi — Yo) + 4V 7? cos? y; sin? (a/2) + 
4V; cos y; sin? (a/2)[AV cos (r — 5 — Yyo)] | 
= V.2+ V2? — 2VoVi cos (yi — yo) + 4V? cos? y; sin? (a@/2) + (as) 


4V; cos y; sin? (a/2)(V» cos y; — Vi; cos ¥;) | 
= V2 + V2? — 2VecVilcos yo cos a cos y; + sin yo sin y;) 
or ABr = By + Bi — 2V BoB; (cos yo cos @ cos y; + Sin yo sin y;) 
Substituting the value of V; from Eq. (6) and setting d(AVr)?/dy; = 0 gives the following equation for y,,, the 
value of y; corresponding to the minimum AV’r: 
2(n — cos Or) tan ym — sin O0e(1 — tan? y,,) 


_ v's — cos 6g + sin Op tan Ym 


= +V B ‘(1 — cos Op) [2 sin yo(n — cos 6g) — cos yo sin Og cos a + sin yo Sin Og tan y,,] (19) 
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Squaring both sides of Eq. (19) yields the following equation in y,,: 


A, tant ym + Ae tan* y,, + A; tan? ym, + Ag tan ym + As = O (20) 
where 
A; sin” Op ‘Ay = [4(y — cos 0g) — D sin? yo sin? Og] sin Op 
A 4(n — cos 0g)” — 2 sin? 0g — D sin yo sin? 0,[5 sin yo(n — cos 0g) — 2 cos yo sin Oz cos a} 
As D sin 6g[cos y% sin Og cos a — 2 sin yo(n — cos Op) |[4 sin yo(r — cos Og) — cos yo sin 6g cos a| — 
4 sin 6e(n — cos Op) 
A; = sin® 06g — D(n — cos 6g) [2 sin yo(n — cos 6g) — cos ~% sin Op cos al? D = By (1 — cos 6p) 


Eq. (20) may be solved for tan y,, in closed form. The method of solution, discovered by Ferrari, may be found 
in any standard text on the theory of equations (see, for example, reference 4, pp. 94-96). As in the previous sec- 
tions, a study of the zeros and infinities of Al’7 and the continuity of the derivatives of Aly with respect to y,; estab- 
lishes that one of the solutions of Eq. (19) gives the minimum Al’y. [See the discussion following Eq. (8). | 

The in-plane minimum (discussed in the previous section) is easily obtained by setting a = (). In general, the 
minimum Al’; (or A8r) is obtained by substituting the y,, obtained into Eq. (18). The corresponding value of 6 is 
found from Eggs. (10) and (11). 

If time of transfer ¢ is important, it might be desirable to use Eq. (18) to plot Al’r versus y; while a plot of transit 


time versus y; may be obtained from the following equations [see Eq. (3); references 2, 5, etc. |: 


@ tan y; 1 — B; cos* y 
sin 06; = —£; cos” ¥; : cos 6; = - 
e e 
VA / cos ¥; e sin (6; + Op) 
= V ab; =e * - 
i ade 2— £6; (1 —e cos (0; + Op) 


. y) ' 
e sin 8, 2 ite. +6 lite. a\Q 
+ a, - | arctan tan — arctan 4 tan (21) 
1 — ecos A, 1 — e l—e 2 l1— 2 
From geometric considerations and from the conservation of angular momentum, the impact variables may be 


found from the following equations: 


—e sin (0; + Op) nB; cos” ¥ 
tan y; = By, = m (22) 
1 — e cos (6; + Op) COs” ¥Y; 

Although the above technique considers only the case of firing when r = r;, the results may readily be modified to 
investigate the merits of firing at a later point in the original orbit. The values of py, eo, and 4) corresponding to the 
given values of 8, yo, and 7; are obtained from Egs. (3) and (21), respectively. However, 6, = 6 + @p where 86, is 
the in-plane range angle from apogee. Therefore, other points along the initial trajectory can be investigated by 
arbitrarily choosing various reference angles 6)’. The corresponding value of #p’ = |)’ — 6,'. Eq. (2) and the law 


of conservation of energy 


E= eo a - = const. 
Zz ri 2 r; 


can then be used to obtain new values of Vj and yo. The lateral range 1 obviously remains constant. The above 
technique can then be utilized for various values of 6)’ and an appropriate firing point can be obtained. 
(E) Transfer From a Circular Orbit to a Point in the Orbital Plane 


An interesting special case of the analysis of the previous section is the transfer from a circular orbit to a point in 
the plane of the circular orbit. For this case, Eqs. (11), (12), (19), and (20), reduce, respectively, to: 


sind = V@,/AB sin ym; cos 8 = V1/AB[VB; cos ¥m — 1] | 
1 — cos Oz 1 — cos Op | 

B= : + 1 — 2 cos ym mere . 
COS ¥m[n COS Ym — COS (Og + Ym)! COS Ym[n COS Ym — COS (Oz + Ym) | 


2(n — cos 6g) tan y, — sin 6p(1 — tan? y») 

y ; =— Vi + cos 6p | 
+ Vn — cos 62 + sin Og tan y,, 

A, tan’ yp, + A» tan® y,, + A; tan? ym, + Aq tan ym + As = 0 


where A, = sin? 0p Ay = 4 sin 6e(n — cos Oz) 
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A;x = 4(n — cos 0g)” — 2 sin? Op A, = —sin Opg/4(q7 — cos Og) + sin Og(1 + cos Oz) | 
; — cos 6 
As = sin? vol me ‘| 
1 — cos Or 


It is noted that y,,, 48 minimum, 6, (/dA/r//”)t, 8;, 8, and y, are all functions of only n and @,. Therefore the 
solution to the problem of the present section can be summarized in graphical form. Figs. 5 through 11 are carpets 
of the above variables versus n and 0g. All the carpets exhibit values for 0g < 180°, and Fig. 6 is extended to values 
of 02 > 180° in order to illustrate that the minimum value of A@,, apparently occurs for 6g = 180°, i.e., for a Hoh- 


mann transfer. 
PART II 


The Error Derivatives 


(A) The In-Plane Problem 
In the previous sections a method was developed for obtaining all V’;, y; combinations for a given total range R 
where 7; and r, are fixed [see Eq. (6)]. Since 7;, r;, and R are specified, and 1’; is a function of y;, the error deriva- 
tives 
OR OR OR Of oO Ot 


, ae ; P , and 
Or; OV; Oy; Or; OV, Or: 
[where ¢ is given by Eq. (21) | can also be obtained as functions of y; alone. 

Let F denote a general function of V;, 7;, 7;, and y; (which may be taken to be either K or ¢). The resulting error 
derivatives OF /O0r;, OF /OV;, and OF/Oy; (r; is treated as a constant throughout Part II) are assumed to be known 
(see reference 6 for the method of evaluating OR/OV’;, Ot/OV;, etc.). 

Since the velocity increment AV; and not the resulting velocity V; is supplied by the vehicle’s propulsive system, 
a transformation of variables from Vj, 7;, y; to AV7, 7;, a, 6 will be performed. This transformation will give the 
error derivatives with respect to the controlled variables. 

Consider first the in-plane case. The transformation from V’,, 7;, y; to AV, 7,;, 6 is given by Eqs. (10), (11), and the 


following equation: 


7 = 7; (24) 
For a given initial orbit, Vy = f(r), vo = a(ri), dVo/dr; = f’, dyo/dr; = gq’. It therefore follows that: 
OAV Vi — Vo cos (yi - 0) 0 =—-«V, sin? (yi — Yo) E ViVi — Vo cos (yi — wt) | 
OV; AV OV;  (AV)?2siné cos 6 (AV)? | 
OF; — OF;_- i OAV V;Vo sin (yi — Yo) 
OV; oo oY, AV | 
06 V/? sin (yi oie I ( ) ViVo sin? (yi — 2] | 
= ~OS i 0 — ~ - 
OY: (AV)? sin 6 cos 6 a bs (AV)? 
OAV _ [Ve -— V; cos (yi — ve) ]f’ — ViVo sin (vi — vo)q’ (25) 
Or; AV 
0 «WV, sin (vi — Yo) [i cos (yi — Yo) — Vo} sin (yi — yo)f" 4 
or; (AV)? sin 6 cos 6 (AV)? 
f ©, Yasin? (ne — wW 7| 
V; cos (yi — Yo) + (AV)? qg 
Of,/Or; = d?,/dr; = 1 
Since 
dAV = (OAV/OV;) dV; + (OAV /Oy;) dy; + (OAV /Or;) dr, 
dé = (05/0V;) dV; + (06/O0y;) dy; + (06/07r;) dr,, etc. 
it follows that 
OF OF OAV OF 06 OF O7, (26) 
= ZU 


Og; OAV dg; 05 Og; OF; dg; 
where g; = Vi, y;or7;. The three equations in (26) can be solved for OF /OAV, OF 06 and OF /OF; to yield the follow- 
ing expressions: 
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OF 06 oAl’ OF 
oF OV, OV; OF OV, OV; 
OAV OOF 8 0 4 86dAV OF 
Oy OY: Ovi OV: po 
D D 
OF OF 7 OF OF OAV OF 06 
OF; Or: ay = const, OF, OAV Or; — 28 Or; | 
6 const 
OF OF O(AV, 4, 7; 
where te and D = : rid) 
Or; Or; V; const ol ] eo Va % ) 


is the Jacobian determinant associated with the transformation in question. 
From Eq. (25) and the elementary properties of determinants (reference 7, pp. 300-306, and reference 8, pp. 22 
26) it is seen that the above determinant D is equal to the Jacobian determinant associated with the transformation 


from V;, y,to AV,6. Its value is given by the following expression: 


OAV 06 
p = (CAV 8 FD) _ O(AT, 8) _ ors OM 
—OVn yer)! |AVnr) ~~ aa Os 
Oy: OY: 
VP aa (vs — We) «2 ,,. Vesin (yi — Ye) 
eR ar — [V;, cos (yi — ve) — Vol= ~~ (28) 
(AV) sin 6 cos 6 AV sin 6 


Evidently [see Eq. (27)] singular points occur when the Jacobian determinant D = 0 or © as is well known from 
transformation theory (see references 9, 10, ete.). From Eq. (28), if AV and V; are assumed finite, singular points 
occur only when 


AV = 0, V; = 0, 7: = Vo or 6=0° or 180 


(B) The General Problem 
A similar analysis will now be presented for the general three-dimensional case where, as before, V) = f(7i), Yo = 


q(r;). From Eq. (15), 


OF OV. = (OF OV;)(0V;/0V,), OF 0a = (OF/OV;)(OV;/0a) (29) 
where OV;/OV, = [2 cos y; sin (a 2) | 
and OV;/0a = —(V, 4) sec y; ctn (a/2) ese (a/2) = (V;/2) ctn (a/2) 
Similarly : 
OF OF OF  £OF OV; OF 4 V, a OF 
= = = esc an y; sec ¥i1=> 
7. Oy vem come. ON OVI OY 2 2 ON, 
a const 
OF Vit OF 
= i te os a 3 
Dy, in 4 or, (30) 


In a manner analagous to that of the previous section, the derivatives OF OV,, OF Oa, OF /O07;, OF Or; are trans- 
formed into OF /OAV7, OF 0&, OF 06, OF O7;, where & = a and 7; = rj. 
It is easily shown that: 


OF 06 OAVr OF 
oF OV, OV, oF OV, OV, 
OAlVr (OF 6.’ 8 DAV, OF 
OV: OF: OF: Ovi (31) 
] J 


OF OF OF OAlr OF 06 OF OF OF OAVr OF 06 


0 Oa AAVr Ia  W2da OF; Oy DAVr 2; 5 dr, 
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OAVr 06 
OV, OV, 
where j= 
OAV, 06 
OV: =O: 


The coefficients can be evaluated with the aid of Eqs. (10), (11), (15) and (18). They are given by the following 


expressions: 


OAV Visec y; — Vc(tan y; sin yo + cos yo cos a) ra 








= ese | 
OV, 2AVr 9 | 
| 
OAVr VV; tan yi — Vo sec y; sin v0) | 
Ovi AVr 
OAV r _ 1 Vie a 3 ; iene e , a = a 
cae iVo ctn 5 cos (yi — Yo) + 2 cos y¥; COS Yo sin’ 5 ¢ — Vctn 5 
AV r yo 7 . sa 
= - ro |; Vo — Vi(cos y, cos yo cos a + sin y; sin ya); f’ + 
Vi Vol(cos y; sin ye cos a — sin y; cos yo)q’] | 
06 — ~‘V,Vo sin’ (yi — Yo) sec Yi iV, —s + ee a ViVe sin? (yi — yo) sec y; a 
- cM LAS > — V; cos (y¥; — ¥e)] csc - = — car 
OV, 2(AV)* sin 6 cos 6 i 2 2(AV)? sin 6 : 2 
06 V2 sin (yi — Yo) | Vo sin (yi — Yo) (32) 
= ——— cos (yi — Yo) + = ——— oe 
Oy: (AV)? sin 6 cos 6 4 %) (AV)? 
) Vo tan y; — Vi sin (yi — ye) — V; tan y; cos (y; — wt | 
06 V;7Vo sin? (yi — Yo) a V7Vo sin? (yi — Yo) a 
= = : V; cos (yi; — yo) — Vo] ctn -- t 
Oa 2(AV)* sin 6 cos 6 \y ” 2 * 2(AV)# sin 6 oe 2 
06 V2 sin (yi — Yo)’ | Vivo ., ( 
= ; - sin 3; = Yo ~— Ce = 0 ae | 
Or; (AV)? sin 6 cos 6 L(AV)? 7 i ia bs 
V,? sin? (vy, — yo)f’ [Ve — Vicos (yi — Yo) ] | 
(AV)! sin 6 cos 6 | 
0a /OV, = 04/07, = 04/dr; = 0; 0a/da = da/da = 1 
O7;/OV, = OF ;/07; = O7;/Oa= 0; O7,/Or, = d¥,/dr, = 1 


The Jacobian determinant associated with the transformation in question is given by the following expression 
(see reference 9, pp. 142-146): 


OA J T 06 Oa OF ; |OAT 'T 06 — 
| OY, OV, OV, OY, | OV, OV, 
OAVr 05 da OF, ‘OAVr 05 a OAV, 8 
as | sll ( a 
OV; OF: OF: OF: Oy: Oy: OV, OV,'| 
=— | = = ( 3 5 
oA Vr 06 0a OF ; OAVr 06 1 0] OAV r 06 | 
| 02 Oa 0a da! | Oa Oa | OF: OF: 
OAVr 06 Ow OF; oA Vr 06 — 
Or; Or; Or; OF; Or; Or; 





i.e., the Jacobian determinant associated with the 


transformation , si ; 
A comparison of Eqs. (31) and (33) shows that the 


1V., Vi a, ri} > {AVr, 4, &, Ft method employed above breaks down at the singular 

points J = Vand J = o. 

i itaaiaaiicaiaicithaes The method of extending the above analysis to ob- 
tain either 0R,/OV,, OL/OV,, ete., [see Eq. (13)] or 


1V., 7:1 — {AVr, 8} the derivatives of R, R,, or L with respect to another 


is equal to the Jacobian determinant associated with 
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Furthermore, it is 
difficult to anticipate which variables will be important 


set of variables is straightforward. 


in a given problem. Any extensions of the above 


theory will therefore be left to the reader. 


FINAL REMARKS 


In Part I a method was developed for obtaining the 
minimum impulse (which results in the minimum fuel 


Velocity and Range Considerations 





Continued from page 684 


(assuming an infinitesimally small attracting body). 
For a range of intercept angles extending to 90 degrees, 
there is a range of range angles which cannot be 
attained. 
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INITIAL FLIGHT TOL 
expenditure) necessary to move a vehicle from an initial 
flight condition to a stationary point target. The 
method requires the solution of a single quartic equa- 
tion. 
investigated by the methods developed in Part I. In 


Other possible modes of transfer may also be 


each case the corresponding values of flight time and 
impact variables (final flight-path angle and velocity) 
are obtained. 

Part II dealt with the evaluation of the error deriva- 
tives corresponding to any chosen transfer orbit. The 
resulting error derivatives are of importance in corre- 
lating initial instrumentation errors with miss distances 
at the target. For a given design problem a trade-off 
between permissible fuel expenditure, flight time, allow- 
able errors, etc., would have to be made with the use of 
the results of both of the above sections. 
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A Theory of Ring-Stiffened Cylinders’ 


R. J. DUFFIN* 
Duke Unwersity 


(1) Introduction 


a ene cylinders are of widespread use in 
present-day engineering practice. For example, 
the fuselage of an airplane is essentially a thin-walled 
cylinder reinforced with uniformly spaced rings. The 
body of a rocket may be characterized in the same way. 
Thus the following problem is of importance in the de- 
sign of such structures. What are the distributed stresses 
when a concentrated load is applied to a ring? This 
paper shows how to obtain ordinary differential equa- 
tions which determine these distributed stresses. 

The first attempt at the solution of the stress problem 
was based on the assumption that the fuselage system 
could be treated as a simple beam. For example, see 
the paper of J. A. Wise.' Experiment showed, how- 
ever, that such a beam analysis must be rather far from 
the truth. A more accurate analysis was given by 
Wignot, Combs, and Ensrud.* They took account of 
the interaction between the ring and the cylindrical 
shell and obtained a differential equation. Unfortu- 
nately, their analysis is not complete. This was pointed 
out by Kuhn, Duberg, and Griffith? in a paper com- 
paring theory and experiment. 

A significant advance was made by N. J. Hoff and 
coworkers at Brooklyn Polytechnic Institute.* ° ° 
The method of Hoff is to express the stress in the cylin- 
der and in the ring as Fourier series in the angle. Then 
the total strain energy is calculated as a function of the 
Fourier coefficients. According to the variational 
principle of Castigliano the equilibrium state is deter- 
mined by minimizing the strain energy. This mini- 
mization process leads Hoff to an infinite set of simul- 
taneous linear equations for the Fourier coefficients. 
Fortunately it proves possible to solve these equations 
by a step-by-step process. The stresses are then found 
by summing a sufficient number of terms of the Fourier 
series. The Hoff analysis is found to give good agree- 
ment with experiment.*® 

In a previous note the writer followed the basic as- 
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terest to note that Griffith does not use variational methods in his 


derivation. 
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sumptions of the Hoff treatment, but Fourier series 
were not employed.’ Instead, classical methods of the 
calculus of variations were used to obtain the Euler 
differential equation. The Euler equation is an eighth 
order ordinary differential equation with constant 
coefficients.** The independent variable is the angle, 
and the dependent variable is the bending moment in 
the ring. The Fourier series method is then seen to 
be merely a special way of solving the Euler differential 
equation. 

An advantage of the differential equation approach 
to the problem is greater flexibility; thus there is no 
difficulty in taking into account a stiffening ring having 
variable cross section. The only difference in this case 
is that one of the coefficients of the differential equation 
becomes variable rather than constant. The solution 
of the differential equation can then be obtained by use 
of digital or analog computers. 

In Hoff’s treatment it is supposed that the applied 
forces lie in the plane of the ring. In a recent report, 
J. Adachi was concerned with the stress problem of the 
ring-stiffened cylinder when the applied forces are 
perpendicular to the plane of the ring.* Adachi ideal- 
ized the problem by supposing that it was sufficient to 
consider an end-stiffened flat sheet. He indicates the 
desirability of having a theory avoiding this idealiza- 
tion. 

In this paper a general approach is developed which 
permits the applied force to have an arbitrary direction. 
Also the assumptions on the properties of the ring are 
Thus in the Hoff treatment only flexure in 
It was 


relaxed. 
the plane of the ring was taken into account. 
assumed that the ring was perfectly flexible perpendicu- 
lar to its plane. In the present theory account is taken 
of flexure perpendicular to the plane and torsion as well 
as flexure in the plane. 

The theory is again based on Castigliano’s principle. 
The strain energy is taken to consist of five separate 
terms resulting from shear stress 7 in the cylinder, nor- 
mal stress o paralle! to the axis of the cylinder, bending 
moment ZL perpendicular to the plane of the ring, 
torsional moment J/ of the ring, and bending moment 
N in the plane of the ring. The Euler equations re- 
sulting may be written as two ordinary differential 
equations of the seventh order relating the moments 


M and N. 


(2) A Cantilever Beam 


In order to be definite, this paper treats the specific 


problem of a cantilever beam. The beam consists of a 
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thin, cylindrical shell of radius a, length ), and thickness 
c. A reinforcing ring is attached to one end of the 
cylinder, and the other end of the cylinder is built into 
a wall, as shown in Fig. 1. The center of the ring is 
taken as the origin of coordinates ¥, J, 2 The = axis 
is taken along the axis of the cylinder. A concentrated 
load is applied to the top of the ring at point A on the 
positive axis. The Z, ¥, and 2 components of the load 
force will be denoted by X), V1, 2. 
load also has a couple, and the #, ¥, and 2 components 
of this couple are denoted by /y, \/,, and N,. 

The angle ¢ measures rotation in the positive direc- 
The initial side of the angle is 
The terminal side of 


The concentrated 


tion about the 2 axis. 
taken to be the positive # axis. 
the angle ¢ intersects the centroidal circle of the ring at 
B. The point B is taken as origin of a rotated coordi- 
nate system x, y, . The x axis is along the terminal 
side of the angle @¢. 

The entire ring will now be thought of as a free body. 
Then the conditions for equilibrium are expressed by 


The z axis is parallel to the 2 axis. 


* . °-. ww . 
the following six equations: 


sa, = +f" 7 sin ¢ dg (la) 
ay, = -f{" T cos ¢ do (1b) 
aZ, = - [sae (Ie) 
lL, = -{" S sin ¢ do (2a) 
M, = -f" S(1 — cos ¢) do 2b) 
N, = -{" 7(1 — cos ¢) do (2c) 


Here S/a? is the force on unit length of the ring acting 
a® is the force on unit 
Of course 


in the z direction. Likewise 7 
length of the ring acting in the y direction. 
S results from a normal stress o, in the cylinder, and 7 
results from a shearing stress 7 in the cylinder. S and 
7 have been defined to have dimensions of moment of 
force in order to simplify the writing of some of the 
equations to follow. 

There will also be a force R/a? on unit length of the 


ting acting in the x direction. Following Hoff we as- 


sume that #& is negligible for thin cylinders, so R is 
omitted in Eqs. (1) and (2). 

In order to avoid unessential complications, the fol- 
lowing three minor restrictions are made: (a) The 
principal axes of inertia of the ring are to coincide with 
the x and z axes. However, the ring does not have to 
have a rectangular cross section as shown in Fig. 1. 
(b) The line of action of the force S is to go through the 
centroidal circle. This will be true if the radius of the 
centroidal circle of the ring is the same as the radius 
(c) The line of action of 7’ is tangent 
This will be sensibly true if 


of the cylinder. 
to the centroidal circle. 
the z thickness of the ring is small compared to its 
radius. 


(3) Stress in the Ring 


To consider flexural and torsional stress in the ring, 
it is useful to regard the portion AB of the ring sub- 
tended by the angle ¢@ as a free body. Such a free 
body is diagrammed in Fig. 2. Acting on the cut at 
B there will be forces Y, Y, Z in the x, y, z directions, 
and there will be couples L, 7, N in the x, y, 2 direc- 
tions. Of course L is the bending moment perpendicu- 
lar to the plane of the ring, ./ is the torsional moment, 
and N is the bending moment in the plane of the ring. 
These quantities are functions of the angle ¢. For 
@ = 0+, we designate these quantities as Yy, Vo, Zp, 
Lo, 40, No. Taking moments at the point B in the 
x, ¥, 2 directions yields the following three equations of 


equilibrium : 
L = Locos ¢ + Mo sin @ — 


> 
aZy sin @ + f Ssin @d&— (3) 
0 


M = My cos @ — Lo sin 6 + 


a) 


aZo/l — cos @) — | S(1 — cos 6) dé (A) 


0 

N = No + aX sin @ + a¥o(1 — cos g) — 
7(1 — cos @) dé (5) 
/0 


Here?@ = @ —E 

It is desirable to write relations (3), (4), and (5) as 
differential equations rather than integral equations. 
To do this the following derivative formulas will be 


employed: 
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? 'p 
- = f S(1 — cos 6) dé, Fo = f S sin 6 dé 
0 0 
? ? 
| f Seos@dé, FF’ = S$ -— f S sin 6 dé 
0 0 


FP’ + F =S 
Here, of course, the superscripts denote differentiation 
with respect to ¢. It is thereby seen that (3), (4), and 
(5) are equivalent to the following three differential 


equations: 
L = —M' (6) 
M'" + M’' = -S (7) 
N'"’ +N = -T (8) 


To apply Castigliano’s principle, it is required to 
have an expression for the strain energy U, of the ring 
as a function of the stresses in the ring. According to 
hypothesis (a) the principal axes of the cross section at 
B are the x and z axes. Let J, and /, be the moments 
of inertia of the cross section at B about the x and z 
axes. Then according to Love (reference 9, § 255, 
§ 258), in the ‘‘ordinary approximate theory”’ the strain 


energy is given by 


U,. = f (LL? + mM? + nN?) dd (9) 
0 


If Eis Young’s modulus and G the shear modulus, then 
/,m, and n are given as 


1 = a/2EI., m= a/2Gh, n = a/2El, 


Here J, is the torsional constant of the cross section 
B. For example, if B is circular, then J, = J, = 
I,/2. 
If the principal axes of the cross section are not 
parallel to the x and z axes, the expression for U’, will 
contain a term of the form LN. This generalization 
will not be treated here. It is also to be noted that 
Eq. (9) assumes that energy due to longitudinal strain 
and transverse shear are negligible. Presumably it 
would be possible to take account of these in the strain 
energy, but it is believed that the ‘ordinary approxi- 
mate theory”’ is sufficiently accurate for most technical 


applications. 


(4) Displacement of the Ring 


The method of this paper is essentially concerned 
with stress and not with strain. However, the theory 
would be somewhat incomplete without relations for 
the strains and the displacements. This section will 
give formulas for the displacements of the ring, as- 
suming that the stresses L, 17, N are known. If only 
the stresses are of concern, this section may be omitted. 

As usual, let u, v, w be the displacements in the v, 
y, 2 directions. Let 8 be the angle of twist in the y 
direction. Then the stress-strain relations corre- 
sponding to the stresses L, 7, and N are 


L/EI, = w"/a? — Ba (10) 
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M/GlI, = B'/a + w’/a? (11) 
N/E, = —u"/a? — u/a? (12) 


To verify any one of these equations, it is noted that if 
the second term on the right vanishes, the equation is 
clearly satisfied. In this case the equation has the 
same form as for a straight beam. 
of the second term on the right can be checked by 
noting that, if the ring is given a rigid displacement, 
the right side of the equation vanishes. 

It is supposed that J and N are known. 
is determined by solving the differential Eq. (12). 
Elimination of 6 between Eqs. (10) and (11) gives 


The correctness 


Then u 


we’ +w’ = —ail"/EI, + a2?M/GI, 


This equation determines w. Then 8 may be obtained 
directly from Eq. (10). The displacement v is deter- 
mined by the equation 


v= —-4 (13) 


This equation merely states that the ring is treated as 
being inextensible. 

The above treatment of the deflection of a ring is 
essentially that given in a paper by H. B. Tabakman.!" 
In that paper Tabakman solves Eqs. (10) and (11) in 
various simple cases of interest. 


(5) Stress in the Cylinder 


Following Love (reference 9, § 343), the so-called 
“extensional solution” of a thin cylinder is determined 
by the equilibrium equations: 


Go, = ap/c (14) 
(Oa,/O0yv) + (O7r/O0z) = 0 (15) 
(Oa,/0z) + (Or/Ov) = 0 (16) 


Here o, and o, are the normal stresses in the y and z 
directions, 7 is the shear stress, and / is the internal 
pressure in the x direction. These equilibrium equa- 
tions may be readily derived by considering a unit area 
of the cylinder as shown in Fig. 3. It is assumed that 
inextensional forces resulting from flexure of the shell 
are inappreciable and can be neglected. 

In the present problem the only external forces on 
the cylinder are at the two ends. In particular the 
pressure pis zero. It follows from Eq. (14) that o, = 0. 
It then follows from Eq. (15) that 7 is independent of 
z. In other words, tr = r(¢). Making use of this last 
relation, Eq. (16) can be integrated to obtain 


o, = —7’s/a + p(¢) (17) 


Here p is a function only of ¢. 
Recalling the definitions of S and 7 and hypotheses 
(b) and (c), it is seen that the following boundary condi- 


tions must hold. 
S = cp/a*, fT = cr/a* (18) 


It was assumed above that there was no radial force 
transmitted from the cylinder to the ring—that is, 
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RING-STIFFENED 


Rk = (0. This assumption is seen to fit in with the ‘“‘ex- 
tensional solution”’ stress state of the cylinder because 
the extensional solution cannot transmit force per- 
pendicular to the surface. 

The strain energy of a unit area of the cylinder is 


given by 
U’, = [(¢,7/2E) + (e,?/2E) + (r?/2G)]e (19) 


Here / is Young’s modulus for the cylinder and G is the 

shear modulus. In the special case treated here, the 

term in o, vanishes. <A discussion of the strain energy 

will be found in Love (reference 9, § 328). 
Substituting the expression for o, given by Eq. (17) 

into relation (19) gives 

clr) cr’ps cp" or 


cl; = = » . 


2a°k ak 2E 2G 


The total energy of the cylinder is given by 


Qn eh 
= fade fv 
J 0 


Carrying out the integration with respect to z gives 


U. = [ h7? + 1(7")? — 271 ‘S + kS? ldo 
(20) 
Here h, 7, 7, and k are constants, defined as 


h = 6/2Ge,1 = 6°/GEca’, j7 = b7/4Eca, k = 6/2Ec 


(6) Displacement of the Cylinder 


As was mentioned above, this paper is primarily 
concerned with stress rather than with strain. For the 
sake of completeness this section treats the strains and 
displacements in the cylinder. If only the stresses 
are of concern, this section may be omitted. However, 
consideration of the strains leads to a better under- 
standing of the problem. 

Now u, v, w will denote displacements of a point on 
the evlinder in the x, y, z directions. Let «, and e¢, de- 
note the normal strains in the y and s direction. Let 


y denote the shear. ‘Then the stress-strain relations 


are: 
Ke, oy, — VO; (21) 
Ee, = o% — V6, (22) 
Gy =T (23) 


Here vis Poisson’s ratio. In the special case considered 


here, o, = 0. The strains are defined in terms of the 
displacements by the following relations (reference 


9, § 334): 


e, = (1/a)(Ov/Od) + (u/a) (24) 
€, = Ow/Oz (25) 
vy = (Ov/0z) + (1/a)(Ow/0¢) 26) 


Eqs. (24), (25), and (26) will now be integrated to ob- 
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tain “, v, and w. From Eggs. (17), (22), and (25) it 


results that 
E(Ow/0z) = —(r’s/a) + p(o) 
Integration of this equation gives 
w= —(r’s*/2aE) + (pz/E) + ¢(¢) (27) 


Here ¢ is a function only of ¢. This relation and Eq. 
(25) are substituted in Eq. (26) and then integrated. 
This yields 

723 p's? ¢'z TZ " 
6a7E iy 2ak “ a r G + 1) vind 


Relation (28) together with relations (21) and (17) are 
substituted in Eq. (24) to obtain the following expres- 


sion for u: 


i : ; 
6a°*k 2ak a G 


, 
vrs avp 
a (29) 


E E 
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If relations (29), (28), and (27) for u, v, and w are 
studied, it will be discovered that all the boundary con- 
ditions at both ends of the cylinder cannot be satisfied. 
For example, Eqs. (13) and (24) are in conflict. This 
difficulty is seen to arise from relations (21) and (22), 
causing e, and ¢, to be proportional rather than inde- 
pendent. To explain such a paradoxical result Love 
introduces the “‘edge-effect”’ (reference 9, § 339). By 
“‘edge-effect”’ is meant that the flexural strains cannot 
be considered as small in the immediate neighborhood 
of the boundary circles of the cylinder. Since the 
flexural strains have been neglected in the above treat- 
ment, it cannot be expected that all boundary condi- 
tions are compatible. 

The expressions derived for u, v, and w are not with- 
out some significance. They do give a measure of the 
displacement in the cylinder away from the boundary 
circles. 

The boundary conditions on stress are assumed to 
hold rigorously. That is, Eq. (18) applies where the 
ring and cylinder meet. If there were another ring 
attached to the other end of the cylinder, a similar 
boundary condition would apply there. 
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(7) The Variational Principle of Castigliano 


The principle of Castigliano may be formulated in 
this way: A linear elastic system 1s in equilibrium under 
the action of prescribed external forces. Then the strain 
energy of the system is a minimum with respect to a vari- 
ation of internal forces. The varied forces must be in 
virtual equilibrium. By “‘virtual equilibrium’’ is meant 
that the varied forces must balance, but the varied 
forces do not have to be compatible with the stress- 
An analysis of this variational 
242. 


strain relationships. 
principle is given in references 6 and 11, p. 
The latter reference may be consulted for techniques 
in the calculus of variations which are now to be 
employed. 

Let U be the total elastic energy of the system. Add- 
ing relations (9) and (20) for lL’, and L’, gives 


if = | [hT? + i(7")? — 2jT'S + RS? + 
/70 


IL? + mM? + nN? \|d@ (30) 


To apply Castigliano’s principle we vary Xo, Yo, Zo, 
lo, Mo, No, S, and 7. Then virtual equilibrium will 
be maintained if we assume relations (1) and (2) for 
external forces and relations (3), (4), and (5) for the 
internal forces. These are subsidiary conditions to be 
imposed on the variation. The variation of l’ may be 


written 


2-16U = [ (aT — i7° + jS'MT + 
/0 


(—j1" + RS)\6bS + IL6L + mMiM + nNiN \do 


(31) 


Here the following integration by parts has been em- 


ployed: 


{ (iT’ — jS)8T'do = 
J0 


2a 


{ (—17"” + jS’)6Tdod + (17" — jS)6T 
0 0 
The second term on the right here is an end-point rela- 
tion. This is recognized as a natural boundary condi- 
tion and may be set equal to zero. 

To analyze the above natural boundary condition 
it is first necessary to note an imposed boundary 
condition resulting from the continuity of shear stress 
7. In particular this means that 7(0) = r(27). It 
is convenient here and in what follows to employ the 
shorthand notation 

Ar = r(2r) — 1r(0) 
Thus the boundary condition imposed by continuity 
may be written as Ar = Vor A7’ = 0. Inspecting the 
natural boundary condition now shows that it may be 
written in the form 


A(i7’ — jS) = 0 (32) 


Consider a variation in which only No, Xo, and Yo are 
varied. Then it results from Eq. (31) that 


SCIENCES 
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0 = f nN|6No + a sin 5 Xp + a(1 — cos @)6Vo ldo 
0 


Since the three variations are independent, the following 
three conditions must hold: 


er 
| nNdo = 0 (:3:3) 
J 
{ nN sin @ dd 0) (34) 
J 0 
| nN cos ¢ dd = 0 (35) 
0 


Similarly, varying Zo, Lo, and J/o yields the sollow- 


ing three conditions: 


350 


[ [—/L sin @¢ + mA/(1 — cos ¢) |dé = 0 
/70 


[ (/L cos ¢ — mM sin dldd = 0 (37) 
/J0 
[ (/L sin @ + mM cos ¢|ldo = 0 (38) 
J 0 


Next we consider a variation when only 7 is varied. 
The value of 6N is determined by the subsidiary condi- 
tion (5). Thus 

clr 
és = -| 67(1 — cos A)dé 
JV 
This expression is substituted in Eq. (31), and the fol- 
lowing identity is employed: 


> 


— / nNbNdé = | nN | 67(1 — cos 0)dddé 
e ( ea ( J 0 


) ) 


| 67 nN(1 — cos @)didd 


Setting the coefficient of 67’ equal to zero gives the fol- 


lowing Euler equation : 
hT — iT” + 4S’ = } nN(1 — cos @#)dé + A, (39) 
Here X, is a function of the form 
A; = a; + b cos $d + Cy Sin @ 


The constants a, 5, and c, are Lagrange multipliers 
resulting from the subsidiary conditions (la), (1b), 
and (2c). These conditions may be written 


f 67° sin ¢ddd = [ 67° cos ¢dé = 
0 0 


f 67 (1 — cos d)dd = 0 
0 


Finally we consider a variation when only 5S is varied. 
The integral of 5/ is transformed in the same way as 
the integral of 6N was transformed above. For 6L the 
following identity is employed: 
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6S sin 0do@ = 


| is | IL sin 6 dé 
J0 Q 


- } ~ IL8Ldt = | IL f P 
0 0 0 


Setting the coefficient of 6S equal to zero in (31) gives 
—jl’ + kS = f [(/L sin @ + 
m.\M(1 — cos 6) |d—é + A» (40) 


Here 2 is a Lagrange function of the same form as Aj, 
but 
(2a), and (2b). 

It is convenient to introduce the following linear 


resulting from the subsidiary conditions (lc), 


combination of relations (39) and (40): 
[" 
0 


Here A, J/, and # are defined as 
K=k—7f/it, J = jh/1, 

The function A; has the same form as ),. 

of the problem is now determined by the integrodiffer- 


KS’ — JT [(/L cos @ + mJ sin @ — 


naN(1 — cos 6)|d—é + A; (41) 
a = jnt 
The solution 


ential Eqs. (40) and (41). 


(8) The Differential Equations 


In general it is more desirable to work with differ- 
ential equations than with integrodifferential equa- 
Also it is desirable to use as few dependent 
To carry out such simplifications 


tions. 
variables as possible. 


let —¢ and —n denote the left sides of Eqs. (40) and 
(41). Thus 

¢ = jl’ —kS (42) 

n = —KS’+ J/1 (43) 


Then successive differentiations of Eq. (40) yield 


© 


¢ | (/L sin @ + mAJ(1 — cos 6@)| dé + Xz 
JV 


%% 
oS 
4 
JV 
oD 
WV 
4 
/0 


Similarly, 
~~ 
=| 
/70 


n 
JV 


|/L cos @ + mJ sin 6 \dé + dq’ 


[—/L sin 86 + mJ cos 0\d—é + Ax” + IL 
successive differentiations of Eq. (41) vield 
(/L cos @ + mM sin @ — 

naN(1 — cos 6) |d—é + XA; 
[—/L sin 6 + mM cos 6 — 

aN sin @ld— + A;’ + IL 
|—/L cos 6 — mJ sin 6 — 


nN cos O6\d—é + v5” + (/L)’ + mM 


Ii 


the above relations \; and A; are of the same form 
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as \;. The following combination of derivatives is 


seen to eliminate the integrals: 


ce’ +e’ = IL)’ +mM (44a) 
n’’’ +9! = (IL)” + (mM)' — aN (45a) 
To achieve greater generality /, m, and m are here 


treated as functions of ¢. Of course they are assumed 
to be differentiable and Making 
Eqs. (42) and (43) yields the differential equations 


periodic. use of 


—j(T +7) + RS? +S) + UL)’ + mM = 0 
(44b) 

K(S® + §™) — J(T® + T™) + 
(/L)” + (mM)'’ — aN ) (45b) 


To reduce the number of variables, S and 7° are 


eliminated from Eqs. (44) and (45) by use of Eqs. (7) 


and (8). This gives the following pair of differential 
equations which determine the solution of the problem: 
(N® + 2N® + N®) — R(M® + 2M + M) - 
IM')’ +mM=0 = (46 

—K(M® + 2M® + M®) + J(N@+ 2N 
N®) — (1M’)” + (mM)’ — aN = 0 (47) 


In this formulation of the problem only the variables 


pertaining to the ring appear. 


(9) The Boundary Conditions 


The basic Eqs. (46) and (47) are of the seventh order 
It is seen that they are equivalent to fourteen first 
order equations; therefore there must be fourteen 
boundary conditions. To take count it is noted that 
Eqs. (1) and (2) give six boundary conditions; rela- 
tionships (33) through (38) yield six more; the 
continuity condition, A7’ = 0), and the natural boundary 
condition, (32), bring the total to The 


first twelve of these conditions are integral conditions 


and 
fourteen 


rather than end-point conditions. 

It is desirable to express all of the boundary condi- 
tions in the conventional end-point form. To carry 
this out, it is seen from Eqs. (3) and (2a) that 


AL. = -f Ssinédé = Ly 
/0 


Making use of Eq. (6) gives 


(48) 


AM’ —L, 


Similar considerations show that the boundary condi- 
tions (2b) and (2c) are transformed into 


AM = AM, 
N, (50) 


(49) 
AN = 


It is seen from Eqs. (5) and (la) that 


la 
f 7 sin & dé = aX, 
0 


AN = 
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Thus Eqs. (la), (1b), and (le) become 


AN” = aX, (51) 

AN” = aY;, (52) 

AM" = az aa M, (53) 

From the integral formula for ¢ and from boundary 
condition (36), it is seen that A¢ = 0. From the 
integral formula for ¢’ and from boundary condition 
(37), it is seen that Ag’ = 0. Similarly, from the 


integral formula for 7 together with boundary condi- 
tions (37), (33), and (35), itis seen that An = 0. From 


Eq. (43) it is seen that AS’ = 0, because An and AT 


both vanish. Eliminating 7’ between Eq. (42) and 


the natural boundary condition (32), it is seen that AS 


= () because Af = 0. Because At’ = 0, it follows that 
A7” = 0. These results are summarized in the five 
relations 


AT = AT’ = AT’ = AS = AS’ =0 (54) 


From the integral relationship for ¢” it follows that 
Ag” = IAL = IL, 


Here boundary conditions (38) and (48) have been 
employed. In a similar fashion the integral relation- 


ships for n’ and n” yield 
Ay’ = IAL = I, 
An” = mAM + A(IL)’ = mAM + IAL’ + l'AL 
An” = (m+))M, — laZ, + IL, 


All fourteen boundary conditions are now expressed 
as end-point conditions. However, for the sake of 
simplicity it is desirable to eliminate the intermediate 
variables ¢ and » and to express the last three condi- 
tions in terms of Sand 7. Thus, since An’ = —KAS’, 


it follows that 
AS” = —IL,/K (55) 
Since AC” = jA7””’ — RAS”, it follows that 
AT’” = (K — kR)IL,/Kj (56) 
Since A7” = 0, it follows that An” = —KAS’’’, hence 
AS’ = [—(m + 1) M, + laZ, — I'I,|/K (57) 


The boundary conditions are now expressed by the 
end-point relations (48) to (57), inclusive. 

For certain degenerate cases to be treated in the 
next section, it is necessary to give two other relations 
resulting from the boundary conditions. From rela- 
tion (44a) it follows that 


Ac’ = —At’ + mAM + IAL’ + l’AL 


Substitution of the boundary expressions already ob- 
tained yields 


Ag’ = (m + 1)M, — laZ, + I'L, 


Since Af = jAT — RAS, it results that 
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AT = [—(m + 1) M, + laZ, — I'L,] X 
(k — K)/Kj (58) 
Likewise it follows from Eq. (45a) that 
An!’ = —An’ + IAL” + 2I'AL’ + 
"AL + mAM’' + m'AM — 7AN 


Substituting the boundary expressions obtained above 


gives 


An!” = (-2/ —m4+/1L, + 
(m + 21')M, = 2I'aZ, _ nN, 


Since An = KAS® + JAT®), it results that 
KAS® = [J(K — k)l/Kj + 272 +m —1”jL, - 
(m + 21’) M, + 2!’aZ, +7N, (59) 


Of course relations (58) and (59) impose no further con- 
ditions on the system. They are merely alternative 


expressions for two of the boundary conditions. 


(10) Four Special Cases 


In many practical situations some of the terms in the 
strain energy will be inappreciable compared with 
others. This will result in reducing the number of vari- 
ables and thereby simplifying the problem. These 
simplified problems could be treated afresh, patterning 
the procedure after the general case just given. It is 
easier, however, to specialize the general solution 
which has already been obtained. This will now be 
done in four special cases of some importance. 

The first case corresponds to the problem treated by 
Hoff. It is assumed that the ring is perfectly flexible 
in a direction perpendicular to its plane. Thus there 
may be no normal stress in the cylinder where it meets 


the ring. So S = 0 characterizes this case. It is seen 
that the variation of 7 again gives Eq. (39). Differ- 


entiating Eq. (39) three times yields the differential 
equation 


—i(T® + T®) + R(T 4+ TO 


)+tuN = 0 


Eliminating 7 by means of relation (8) gives the fol- 
lowing differential equation for the bending moment in 
the ring 


i(N® + 2N® + N@®) — 
k(N® + 2N® + N@) + nN =0 (60) 


Three boundary conditions are given by Eqs. (50), 
(51), and (52). Thus 


AN = MN, AN’ =aX,, AN” =<ayYF, (61) 


Only forces perpendicular to the axis of a cylinder are 
considered. Similarly, only couples are applied parallel 
to the axis of the cylinder. Thus Z; = L; = M, = 0. 
The boundary conditions are obtained from Eqs. (54), 
(56), and (58) ; thus 


AT = AT’ = AT” = AT® = AT® =0 (62) 
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Together, Eqs. (61) and (62) account for eight bound- 
ary conditions for the differential equation (60) of the 
eighth order. 

The second case goes to the opposite extreme. Thus 
it is now supposed that the ring is rigid perpendicular 
More precisely this case is characterized 
This means that there is no 


to its plane. 
by taking / = m = 0. 
strain energy resulting from flexure in the z direction 
or torsion in the y direction. Let y = K¢’ — kn = 
iKT” — RJT. 

It follows from Eqs. (44a) and (45a) that y’/"" + y’ = 
knN. Then 7 is eliminated from this relation by 
means of Eq. (8), resulting in the eighth order differ- 
ential equation 
—jK(N® + 2N® + NM) + 

kRJ(N® + 2N® + N@) — kaN =0 (63) 
The boundary conditions are the same as before—that 


However, relation (62) now follows 
It is not necessary to assume 


is, (61) and (62). 
because / and m are zero. 
that Z,, L;, and Wf, are zero. 
The third case to be considered concerns a ring which 
is flexible in its own plane. The applied force is con- 
This case is char- 
Thus Eq. (46) be- 


sidered to be parallel to the z axis. 
acterized by the relation N = 0. 
comes 

k(M® + 2M + M@) + (1M')’ —mM=0 __ (64) 
Boundary conditions (48), (49), and (53) apply. Thus 


AM = MN, AM’ = —Ly, AM” = az = MN, (65) 


Boundary conditions (54) and (55) give 


AS =0, AS’ =0, AS” = —/L,/K (66) 
Thus Eqs. (65) and (66) are six boundary conditions 
for the sixth order differential equation on the torsional 


moment. 

The fourth case is the opposite extreme to the third 
case. Thus it is now supposed that the ring is rigid 
in its plane. This case can be characterized by the 
relationn = 0. Let 


QO = Jt — jn’ = jKS" — kJS 
Making use of (44a) and (45a) gives 
O'’ + Q' = —j(IL)'"’ — j(mM)" + JIL)’ + JmM 


Eliminating S by means of relation (7) yields the fol- 
lowing eighth order differential equation on the tor- 
sional moment 


—jK(M® + 2M® + M®) + kJ(M® + 
oM® 4+ M®) = j(IM’)!" — j(mM)" — 
J(1M’)’ + JmM_ (67) 


The boundary conditions include (65) and (66) and the 
two further conditions (57) and (59). 


(11) Further Generalizations 


(A) The stress distribution resulting from a concen- 
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trated load applied at the angle ¢» is, of course, obtained 
by rotating the standard solution just found through 
the angle q@. 
corresponding to a system of concentrated loads is ob- 


More generally the stress distribution 


tained by superimposing such rotated solutions. 

(B) The differential equations corresponding to a dis- 
tributed load are given by similar methods. More- 
over, the standard solution obtained here serves as a 
Green’s function for the continuous loading problem. 

(C) The differential equations for a cylinder with 
several reinforcing rings are obtained in essentially the 
same way. It is seen that a system of equations re- 
sults involving three moments L,, .1/;, and N; for each 
ring and two stress variables S; and 7; for each section 
of the eylinder. 

(D) The cylindrical shell does not need to be iso- 
tropic. Hoff made use of this notion to incorporate 
the effect of stringers—that is, longitudinal stiffening 
rods. Likewise the treatment applies if the cylinder 
is corrugated, as in the problem studied by Adachi.* 

(E) Itis permissible to have the parameters vary with 
the angle ¢. 
without change, for a ring of variable cross section 


The above treatment applies directly, 


because it has not been assumed that the coefficients 
/, m, and m are constant. Moreover, the treatment 
could be modified to take into account a cylinder 
whose thickness varies with the angle ¢. 

(F) A similar theory could be developed for cylinders 
which are not circular, such as elliptic cylinders. 
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Vortex Containment for the Gaseous-Fission 
Rocket’ 


JACK L. KERREBROCK* ann ROBERT V. MEGHREBLIAN** 
California Institute of Technology 


Summary 


The nuclear rocket is potentially capable of much higher 
specific impulse than any chemically fueled rocket, because of the 
high energy content per unit mass of the fissionable material. 
While a part of this potential can be realized by use of a low 
molecular weight propellant heated by solid fuel plates, it seems 
clear that the full potential can be realized only if the fissionable 
material can be retained in gaseous form, and its fission energy 
transferred directly to the propellant. 

The purpose of this paper is to consider the possibility of re- 
taimug the fissionable material, in gaseous form, in the pressure 
field of a vortex, while the low molecular weight propellant dif- 
fuses radially inward through it. A two-dimensional, laminar, 
vortex flow is studied, the propellant being introduced tangen 
tially at the periphery and withdrawn at the center. 

The mass-flow capacity of the vortex, per unit of vortex length, 
is found to be independent of the vortex diameter, so that a 
large number of small-diameter vortices filling a given volume 
will have a much larger mass-flow capacity than a single vortex 
of the same volume. At tangential Mach numbers of the order 
of unity in a mixture of hydrogen and plutonium, the mass-flow 
capacity is of the order of 0.01 pounds per second-foot. 

If the vortex pressure gradient decreases monotonically out 
ward, it is possible to maintain the heavy gas in a highly peaked 
distribution, away from both the center and the periphery of the 
vortex; however, for very high heating rates, or very large tem- 
perature ratios, the pressure gradient cannot decrease mono 
tonically outward, and in this case the concentration of heavy 
gas near the periphery of the vortex will be high. 

Insofar as the vortex-containment mechanism limits the gase 
ous nuclear rocket, its performance potential is very high. It 
may be that its performance will actually be limited by other 
factors, such as radiant heat transfer from the gas to solids, the 
difficulty of maintaining fissionable material in gaseous form, 
and the difficulty of generating vortices with the required low 
radial mass-flow rates and high tangential velocities. The po- 
tential of the gaseous fission rocket seems to warrant intensive 


study of these problems. 


Introduction 


i lew NUCLEAR ROCKET, in which the fission reaction 
acts as a source of energy to heat the propellant, 
which is then expanded through a nozzle to produce 
thrust, is potentially capable of much higher specific 
impulse than is any chemically fueled rocket. This 
is because the energy content of the heated propellant 
per unit mass is not jimited, as in the case of the chem1- 
cal rocket, by the chemical energy content of the pro- 
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pellants. The energy content of the fissionable mate- 
rial per unit mass is in fact so high that it need not be 
used as propellant; its energy may be extracted and used 
to heat a separate propellant. 

However, the problem of transferring the energy 
from the fissioning fuel to the propellant makes full 
realization of this great potential difficult. A portion 
of it can be realized by choosing a propellant with a 
very low molecular weight, so that even with the pro- 
pellant temperature limited to the values obtainable 
with heat transfer from fissionable materials in solid 
form, the propellant enthalpy, and hence the specific 
impulse, is considerably above that attainable with 
chemical propellants. It seems likely that, with hy- 
drogen as the propellant, values of specific impulse 
ranging from S00 seconds to perhaps 1,200 seconds 
may be obtained, depending on the pressure at which 
the hydrogen is heated. 

It seems clear, though, that the full potential of the 
nuclear rocket can be realized only if the fissionable 
material is retained in gaseous form, and heat is trans- 
ferred from it to the propellant without the interposi- 
tion of any solid surface. Two methods are available 
for this energy transfer: radiation and direct molecular 
collisions. Of these, we shall consider only the latter. 
To effect the direct transfer of the kinetic energy of the 
fission fragments to the propellant by molecular colli- 
sions, the fissionable material and propellant must be 
in intimate gaseous mixture within a nuclear reactor 
in which a chain reaction can be supported. It is 
readily demonstrated that the requirement for nuclear 
criticality, together with the desire for large amounts 
of heat addition to the gaseous propellant, dictate 
rather high concentrations of fissionable material in the 
gas mixture. If this mixture is then simply exhausted 
through a nozzle, the loss of fissionable material is 
prohibitive from the standpoint of cost. Thus, it is 
essential that the fissionable material be retained within 
the reactor, while the propellant flows through it. 

The purpose of this paper is to consider the possi- 
bility of containing the fissionable material within a 
vortex flow. The strong radial pressure gradients 
which can be generated in vortices are capable of pro- 
ducing appreciable diffusion of a light propellant gas, 
such as hydrogen, through a heavy gas, such as plu- 
Thus, it seems possible that a distribution of 
vortex while the 
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mentation of this concept is shown schematically in 
Fig. 1. The propellant gas, and enough fissionable 
gas to make up losses, would be introduced with a high 
tangential velocity at the periphery of the vortex tube, 
through a series of jets distributed over its surface. 
The propellant would then diffuse radially inward, 
through the distribution of fissionable gas, and leave 
the tube axially near its center. 

For the purposes of analysis, this flow will be idealized 
to a two-dimensional vortex flow, into which the fluid 
is introduced uniformly at the periphery with a tangen- 
tial velocity, and from which it is withdrawn at some 
radius near the center, as by a cylindrical sink at that 
radius. No attempt will be made to treat the three- 
dimensional flow near the axis; it is felt that this flow 
is not of crucial importance to the containment mecha- 
nism. A question does arise as to the proper boundary 
conditions to be applied at the inner and outer radii 
of the two-dimensional flow; this question will be dis- 
cussed in some detail. 

The flow will be assumed to be laminar. 
that this assumption is rather extreme, because the 
Reynolds number of the flow, based on the high tan- 
gential velocity and the tube diameter, must be large, 
and the resulting turbulent mixing would tend to in- 
validate the laminar calculation. However, the na- 
ture of turbulent vortex flows is not well understood, 
particularly when the fluid is compressible, so that 
conclusions based on turbulent mixing processes for 
radically different flows (e.g., pipe flow) may be com- 
pletely erroneous. The point of view which is taken 
here is that the validity of the laminar diffusion calcu- 
lation can best be judged by direct comparison with 


It may seem 


experiment. 

Fission-fragment heating of the propellant is as- 
sumed to occur locally at the point of fission—1.e., 
the range of the fission fragments is assumed small 
compared to the radius of the vortex. At one hundred 
atmospheres pressure and 2,000°K., the range of fission 
fragments in hydrogen is about 0.7 cm.,? however, half 
of their energy is deposited in the first 0.2 cm., hence, 
at this pressure and temperature, the vortex radius 


COOL GAS AND FISSIONABLE MATERIAL 
INTRODUCED TANGENTIALLY 


must be one centimeter or larger. An additional as- 
sumption, that the heating rate is proportional to the 
concentration of fissionable material, implies that the 
neutron flux within the vortex is constant. It will be 
assumed that the mole fraction of fissionable gas in the 
mixture is small. This assumption simplifies the anal- 
ysis considerably, and because of the much larger 
mass of the fissionable atoms compared to that of the 
propellant atoms, is valid even when the density of 
fissionable material is considerably greater than that 
of the propellant. 

Finally, the vortex flow will be assumed inviscid, so 
that the tangential velocity is inversely proportional 
to the radius. The justification for this is that the 
radial-inflow Reynolds number is sufficiently large for 
the interesting cases that viscous effects are in fact 
small, except near the solid outer boundaries of the 
vortex. Boundary-layer-like flows with high rates 
of shear will exist on such surfaces, and the resultant 
torque may make the attainment of vortices with 
sufficiently high tangential velocities difficult ; however, 
this problem will not be considered here. 

In the present paper, we shall be concerned with the 
process of diffusion of the propellant through the heavy 
gas, with heat addition. No detailed consideration 
will be given to the many other problems, such as 
radiant heat transfer and nuclear criticality; which 
must clearly be of great importance for a gaseous 
reactor. 

Nevertheless, some estimates of the performance of 
gaseous nuclear rockets, as limited by the vortex con- 
tainment, will be given. It must be recognized that 
these estimates cannot reflect the ultimate potential 
of the gaseous nuclear rocket with vortex contain- 
ment. They can only serve to indicate the limits 
imposed by the containment scheme itself. 


Formulation of the Diffusion Problem 


For a complete description of the two-component 
gas mixture, two new variables must be added to the 


usual fluid dynamic set. These may be taken as the 
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Fic. 1. Specialization of the vortex tube for rocket propulsion, 
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density and velocity of one species, say the heavy gas, 
in which case the continuity equation for the heavy 
gas and the relation connecting its diffusion velocity 
to the pressure and concentration gradients complete 
the set of governing equations. 

In the two-dimensional vortex flow, the continuity 
equation for the mixture is quite simple. If u is the 
radial gas velocity, p the density, and r the radius, then 

pur = (IM, + Me) /2r (1) 
where 9M, and IW» are the constant mass flows of the 
light and heavy gases, respectively, per unit of vortex 
length. Similarly, since the two gases do not react 
chemically, the heavy gas satisfies the continuity 
equation 


pu + us)r = Me/2r (2) 


where 7 is the diffusion velocity of the heavy gas. 
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Fic. 2. Variation of density ratio, «, with dimensionless radius, 
r*, for g,, less than g,, (crit) 
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The Diffusion Equation 
Neglecting thermal diffusion, the diffusion velocity 
is related to the concentration and pressure gradients 
by :? 
p f n_ d log (n2/n) 
U» = —D»- — 
p — pe a nN» dr 
n(m, — mz) d log p\ (3) 
as. . 
p dr J 
where 7 is the total gas concentration, p the total] gas 
pressure, m, and my, the respective molar masses, and 
Dy, is the binary diffusion coefficient. Now by using 
Eq. (2) to eliminate #2, Eq. (3) may be put in the form 


( Me s *) —2arp Dy (“) x 
M, + Me p MW, + We \p 


f n_ d log (n/n) n(m, — m2) dlog p\ 
p dr 


(4) 
— N» dr 
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It is convenient to rewrite this expression in terms of 
the ratio of fissionable material and propellant densi- 
ties, w = po/p;. As was mentioned previously, it 
will be assumed that m./n < 1, however, because m, >> 
m,, it does not follow that w < 1. In terms of w, 
the diffusion equation then becomes 


( Meo Ww ) —2rrpDy ( Ww ) 
MM, + Me l + Ww MM, + MWe ] + Ww 


fdlogw = (m2/m, — 1) d log p| =) 
(dr l -- ww dr § 


It will become clear later that I», 9, is determined 
by the value of w at the boundary where the flow leaves 
the vortex (in this case, the inner boundary). Thus, 
M2/IM, is a characteristic value of the diffusion equa- 
tion. The magnitude of this characteristic value in- 
fluences the behavior of w where it is of the order of 
MW./IN,; where w > Me/My, the behavior is controlled 
by the factor 2mrp,Dy2/(9, + M2), and by the be- 
havior of the pressure gradient. 

The radial flow velocity must be of the order of the 
diffusion velocity, and it is readily verified that it is 
therefore small compared to the tangential velocity. 
In this case, the radial momentum equation becomes 
simply 

dp/dr = pv*/r 


where v is the tangential velocity; using the fact that 
nN» << n,, this may be written: 


(d log p)/(d log r) = (1 + w)v?/RT (6) 


where F& is the gas constant for the light gas. Substi- 
tuting this result into Eq. (5) and multiplying through 
by (1 + w)/w yields 


a (' + “\( Me )- 
w /\O +m 


2rp,D2 k log w (™ = ) vy? | (7) 
IM, + Me Ld log r my RT 
The most important features of the vortex contain- 
ment mechanism can be seen quite readily from Eq. 
(7). Suppose first that there is no radial flow, so that 
MM, = Wy. = OV; then the quantity in square brackets 
must be zero, and we have 


l (™ ) I ( : )¢ +, tant 
o¢g WwW = _ % Ste 
gu “| RT) > constan 


Thus, w increases monotonically with r for mz/m,> 1. 
This is of course the equilibrium distribution for a given 
variation of v and 7 with r. 

Next, suppose that 31, and IN, each have some nega- 
tive value, but that || << [Ny |: then, except where 
w is very small, the left side of Eq. (7) is unity, and 
except for the variation of the quantity 27p,D.» 
(MN, + Mz), the sum of the two terms in brackets must 
be a negative constant. Thus, initially as 31, decreases 
from zero, dw/dr must decrease. If 9M, is made suffi- 
ciently negative, dw/dr can be made zero at any given 
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value of r. Now, for such a value of 3%, if v2’ RT de- 
creases monotonically with increasing 7, dw/dr will be 
positive for 7 less than the point at which dw/dr is zero, 
and negative for y larger than this value. Thus, a 
simple maximum is produced. 

The position of this maximum may be obtained by 
solving Eq. (7) with dw dr = 0. Ti Mim, = v/(yRT)'” 
is the tangential Mach number at the point of the maxi- 
mum, and other quantities evaluated at the maximum 
are similarly denoted by subscript m, then 


: coms ] NM, + Me 
y¥M im? = x 
(m2/m,) — 1 2rpDy /m 


l Wm We 
[ ~ ( = Y( | (8) 
Wm MM, + My 


Eq. (S) contains one of the most important pieces 
of information to be derived from the analysis, namely, 
that the radial mass flow permitted by the diffusion 
process is independent of the diameter of the vortex, 
and in fact depends only on II, vy, AZz,? and the 
product (p:Dj2),, evaluated at the maximum point of 
w, since for the cases of interest SW2/(, + WM) <« 
W, (1 + W,). For a given pair of gases, (p,D12) 
depends only on 7°,,.. The most important fact is that 
the mass flow rate per unit of vortex length is inde- 
pendent of the vortex radius. Physically, this is a re- 
sult of the fact that for a fixed Mach number, the pres- 
sure gradient is inversely proportional to the vortex 
radius, and the resulting diffusion velocity is therefore 
also inversely proportional to the radius of the vortex. 
But the mass flow rate per unit of vortex length is pro- 
portional to the product of the radial velocity (which 
is essentially the diffusion velocity) and the vortex 
radius, and hence is constant 

Finally, if v?/ RT is not monotonic but first decreases 
and then increases with increasing 7, w can have both 
a maximum and a minimum. It will be demonstrated 
that this situation can occur with heat addition in a 
free vortex. There will then be two positions in the 
vortex where dw/dr = (. In a free vortex, where 
v«1/r, there will be a level of heat addition below which 
only a simple maximum will exist because 7 does not 
vary rapidly enough to overcome the variation of v. 
This value will be termed the critical rate of heat addi- 
tion. There will also be a larger rate of heat addition 
above which no maximum of w can exist, because 
(v2 RT) increases monotonically with increasing r. 
This value will be termed the maximum rate of heat 
addition. Both of these values are of course de- 
pendent on 27p,Dy2/( + Me), Me/IMy, m2/m, and 
M im. 


The Energy Equation 

To make these relationships quantitative, we must 
consider the energy equation, since it describes the 
behavior of 7. The viscous dissipation will be neg- 
lected, as it is quite negligible in a flow with large heat 
addition. The equation in terms of the enthalpy h is 


then® 
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dh dp ld 
pu =U _ 


(7g) + 2 (9) 
dr dr r dr t . 


where g = A(dT/dr) + (Cr — Cy) pet27, X being the 
heat conductivity, and g is the volumetric rate of heat 
addition. In this case, g of course represents the heat- 
ing by fission fragments. 

The heating rate is proportional to the concentra- 
tion, ms, of fissionable material, to its fission cross 
section, o,, and to the neutron flux, ¢. If the local 
energy release per fission is Q, then g Qo; on. In 
the present work, it will be assumed that o, and ¢ are 
independent of 7. Q is a constant for practical pur- 
poses, so the heating rate is simply proportional to 7». 
h is the enthalpy of the gas mixture, and under the 
assumption that m). << n, it can be written ash = C,,T 
(1 + w). The specific heat of the light gas will be 
taken constant, although this is a very crude assump- 
tion for diatomic gases at the temperatures which are 





12 ~ +————_————_y —y 
Mtm=0.7) Wm=l.0 a 
| 9m /Gm (MAX) 
OF | 
0.8} | 
T° | 
06} 
0 ; } 
= 
ee ee ee 2 Ee ee a hued 
0.6 07 0.8 09 1.0 Li 1.2 3 


Fic. 3. Variation of dimensionless temperature, 7*, with dimen- 
sionless radius, r*, for g,, less than g,, (crit 


of interest for the gaseous fission rocket. We make 
the assumption in order that the calculations shall not 
depend explicitly on the temperature and pressure 
level, as they would if the exact relationship between 
h and T were included. With the expression previ- 
ously given for dp/dr, Eq. (9) then becomes 


d ( CnT ) v? 24 
= - meee 
dr l a Ww r IM, a Me 


1 dT 
- g “ al (Cre —_ Cures | + Qa; QNeo pu 


dr dr 
From Eq. (2), 
polor = 
[(My + Me) /2r] [Me/(M, + Me) + w/(1 + w)) 
and from Eq. (1), 
Qo; on2/pu = [27Qo; o/m2(M, + Me) |rpw 


so that 
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d : v" Cpe related by \ = (9y — 5)uC, 4, so that the ratio becomes 
oe +{1— x 


9y — 5 (’ 4 “( 27 ) (11 
d ( Mee es ) r| be 2a x 4y w My + Me 
IN, + Me I + Ww MM, + Me 
d ( \ =) 2nrQa; o 
r + : : 
dr Cc moCy (IN + M2) 


It will become clear later that w must be of the order 
me C0) of unity or greater over most of the vortex to provide 
dr sufficient fissionable material for nuclear criticality, 
so that the factor (1 + w)/w is of order unity. The 
factor (9%, + My) 2ayu is the radial Reynolds number 


of the flow; for the situations of interest it is of the 


pl 


It will be helpful to consider the relative magnitudes 
of the heat transfer by conduction and diffusion. As- 
suming that rd7‘/dr and 7 are of the same order, the 
ratio of the conductive and diffusive heat transfer rates 
is [27A/C,, (9%, + Me)|/[w/(1 + w)], at points where 
w >> Me/M, Now according to Euchen’s formula,’ 
the heat conductivity, viscosity, and specific heat are 


order of 100. Thus, conduction can be neglected com 
pared to diffusive energy transport, except where 
r(dT/dr) > T. 

Neglecting the conduction term, carrying out the 





indicated differentiation, and eliminating dw/dr with 
Eq. (7), we find the following expression for (d log 7) + 


(d log r): 
l+w 2rOo; or*piw 
: c — 1) a ll 


d log T - w mC, (MN, + M2)T 4 


d log r l Cre is ee Me 
+(1-"}ii- 
w Can w WM, + Me 
Cy l My + Me l+w Me Ms 
gree) anda L'~ (oat a) + Gos — 1) 4 
Ca \1l + 'w ] 2rp\ Do w MM, + Me my f 
l Cos l W Wie 
+ (1 EN OS") a) 
w Ch w IM, + My 
From Eq. (7), the accompanying expression for (d log w/d log r) is 


log w | Me l w Me » 
d log 7 ( i+ aE _ ( + )( ly )| ~ (™ _ ') Mt? on 
d log r 2a pyDy» w MM, + Me my 


(d log p)/(d log r) = yl/7(1 + w) (14) 


and from Eq. (6), 


With the proper boundary conditions, and the relation p; = p 7’, these three first-order nonlinear equations may be 
solved for 7, w, and p, if v is given as a function of r. As was mentioned in the introduction, it will be assumed in 
the present work that v « 1 /r, i.e., that the flow is inviscid. 

Only those solutions yielding a maximum of w within the vortex will be considered. For heating rates below the 
critical heating rate, such solutions offer the possibility of containing the fissionable material within the vortex, away 
from the outer wall. For heating rates above the critical value, the value of w at the outer boundary of the vortex 
will be higher, in contrast. 

To facilitate the numerical integration, it is desirable to rewrite the equations in terms of dimensionless variables. 
A particularly convenient set of variables is obtained by dividing each of the variables by its value at the radius, 
’», Where wis a maximum. Such variables will be denoted by an asterisk—e.g., 7* = r r,, T* = T/T». It is also 
convenient to introduce a dimensionless heating rate g,,, evaluated at 7,,, as 


m= = [21Qo,; PP in" Nom ] [Cn ( My = Jz) Ts ( I 5) 


Through Eq. (8), (SM, + M2) /2rp; D2 may be related to AZ;,, w,,, and My (IM, + Mz), provided the dependence of 
D,» on the state of the gas mixture is specified. For the present purposes, the simple result for Dj. based on hard- 
sphere collisions, namely, 


Dy = (3/S8ndy2") [(m, + me)RT /2armymz)!/? (16) 
isadequate. kis Boltzmann’s constant and d); is the effective collision diameter. Thus, Dj, (Dij2)m = T**/? p*, and 


) 


mM, + mM. _ (= tf = | Tein YM m?((m2'm,) — 1|T*- 


2rpDy» 2rpD 2 he (' re =) Me (17) 
Wm M, + Me 


With these relations. Eq. (12) becomes 
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l+w | — 1)M,,,” ead 
— gm 


d log ‘fag Ww pr 1 **w,, 


dlogr* 1 +(1 _ eV _ (i+ “\( IM, 
w ee w M; + Me 
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_ 


r (= *\ My» ) 
YM tm?[(m2/m,) — 1](Cpe Cpr) l Ww Mm, + My 


(1 + w)7T*!/: 


Similarly, Eq. (13) 1s 


YM tm?[(m2/m,) — 1} 


d log w 


d log r* 


and Eq. (14) reduces to 


T*} 2 r*T*! 2 





p*2 7 *1/2 si (‘ + “=)( We ) 
= Wm NM, + Me 
. ‘e ~ ( 1S) 
+ (; oo -) i (’ + “\( Me ) 
w Cask. w Mm, + mM 
' (‘ + “\( M:. ) 
w MM, + Me 
(19) 
' (| Me ) 
Wy, MM, + Me 
YM im?(1 + w)/r*?T* (20) 


(d log p*) (d log r*) 


Boundary Conditions 


The system of three first-order differential equations, 
(18), (19), and (20), allows one boundary condition 
to be satisfied on each of the three variables w, 7*, and 
p*. In addition, the quantity IN./(I, + My) is a 
characteristic value of the diffusion equation, as will 
be demonstrated. It is of course clear that by neglect- 
ing conduction in the energy equation Eq. (18) is re- 
duced from a second-order equation to a first-order 
equation; this reduces the number of applicable bound- 
ary conditions to one. In a vortex such as is sketched 
in Fig. 1, the boundary condition on 7 at the vortex 
tube surface, namely, that the gas temperature equals 
the surface temperature, is lost. Similarly, by neglect- 
ing the radial acceleration in the radial momentum 
equation, Eq. (20) is reduced to a first-order equation. 
If the pressure is specified at any radius in a tube such 
that in Fig. 1, Eq. (20) appears to determine the pres- 
sure everywhere else. Actually, the three-dimensional 
flow near the axis of such a tube would be an important 
factor in determining the pressure level in the tube. 

The situation is perhaps not so clear in the case of 
the diffusion equation. If, for example, the vortex 
is assumed to be bounded at its inner and outer radii 
by porous cylinders, it may appear that the value of 
the density ratio, w, in the vortex just adjacent to the 
surface of each of the porous cylinders should equal 
the value of w in the gas mixture flowing through the 
porous cylinder. However, this is only true if the 
direction of gas flow is from the vortex to the porous 
surface, i.e., if the boundary is an outflow boundary. 
At such a boundary, the mass flow rate of each species 
is proportional to its concentration. The inner bound- 
ary of the present vortex flow is such a boundary; if 
the radius of this inner boundary is 7,, and the corre- 
sponding value of w is w,, then it is clear that the quan- 


tity WM. (MM, + M.), which is the ratio of heavy gas 


mass flow to total mass flow, is simply 


Me (MN, 2 We) = w,. (1 + W,.) (21) 


Thus, the characteristic value of the diffusion equa- 
tion is determined at the inner boundary. 

In contrast, all that is required at the outer, or in- 
flow, boundary is that the gas flowing into the vortex 
have the same w as that leaving, namely, w,. It is 
neither necessary nor in general possible for the value 
of win the vortex adjacent to the inflow boundary to 
equal w,. The transition from one value of w to the 
other at the inflow boundary occurs by a mixing proc- 
ess, and results in an increase of the entropy of the gas 
mixture. This entropy increase due to mixing is 
offset by the separative effect in the vortex, so that 
the entropies of mixing of the inlet and outlet flows are 


equal. 
For the present calculations, the boundary condi- 
tions will be applied at r* = 1, and are simply 


7*(1) = 1, 


p*(1) = 1, wl) = w, (22) 


where w,, is the maximum value of w. For given values 
Of Lm, We, Wm, Mim, M2/m,, Cyo/ Cp, and y, the solution 
can then be obtained by forward integration. The 
inner radius of the vortex is then defined as that radius 
at which w is equal to w,. 

The last statement implies that the fluid leaves the 
vortex at the radius 7,, as though there were a cylin- 
drical sink there. This condition would be approxi- 
mated by a porous cylinder of radius r,. The three- 
dimensional flow near the axis of a tube such as that 
sketched in Fig. 1 would act as a sink distributed over 
the portion of the tube volume within the radius of the 
exit nozzle. For a precise evaluation of the loss of 
fissionable material from the vortex, a treatment of this 
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Fic. 4. Variation of density ratio, w, with dimensionless radius, 
r*, for g» greater than g,, (crit). 


distributed sink would be essential; for the present 
purpose, w, will be taken to be very small (10~*), and 
the fluid will be assumed to leave the vortex at 7. 
This gives a very conservative estimate of the loss 
rate. 

Finally, it should be noted that for w,/w,, << 1, the 
* 


character of the solutions for 7*, p*, and w is inde- 


pendent of the magnitude of w, [see Eq. (19) ]. 


The Critical and Maximum Heating Rates 

The maximum heating rate has been defined as the 
highest heating rate for which w has a maximum within 
the vortex. Equivalently, 2n (mar) is the largest value 
of g,, for which w has a maximum at r* = 1; that is, 
the largest value of g,, for which d°’w dr*? is negative 
at r* = 1. From Eq. (19), 


d*w My 
_ = M ae — — l x 
(=). _— (“: ) 


Wm [2 + (1/2)(dT*/dr*),+=1] 
so that for 2mimar), (d7 */dr*),*«-,; = —4. Now from 
Eq. (18), 


l + Wm [(y —_— 1) Min? eins £m] 


(=) Wm (23) 
eat>) 

ir* | ,¢= l Cre 1 — (W./Wm 

6 l + (; = I ) ( ) 


Wm 
whence 


2m(mar) = (Y uae 1) M tm? + 


l Cro\f 1 — (We/Wm Wy 
4 | ~ (: — — Y( | . ) (24) 
Wm C pl l + W, 1 + wv, 


and forw, « land w, K w,,, this is 
Zm(mar) = + = (y = 1) M tm? 


The critical heating rate has been defined as the 
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smallest heating rate for which dw dr* = 0 for some 
r*> 1. From Eq. (19), for w. K wy, dw/dr* is zero 
if 7* = 1 r**; however, since 7* can be obtained only 


from integration of the differential equations, the value 
Of Lm(crit) Cannot readily be obtained analytically. The 
values Of 2m(crit)’ Zmimar) Obtained from the solutions of 
the equations will be given later. 


Solution of the Diffusion Problem 


Solutions of a system of equations equivalent to 
Eqs. (18), (19), and (20) with boundary conditions 
(22) were given in reference 1. These were obtained 
by forward integration of the equations, from r* = | 
to both smaller and larger 7*. The parameters which 
characterize the gas mixture, and which occur ex- 
plicitly in Eqs. (18), (19), and (20), were chosen to be 
representative of hydrogen and plutonium gas. Spe 
cifically, m2/m, = 119.5, Cp/Cy = 0.008, and 
1.31 were chosen. The exit density ratio, w,, was 
taken as 0.0001. 

The remaining parameters are then g,,, which meas- 
ures the heating rate; .\/,,,, which measures the strength 
of the vortex; and w,,, which is a measure of the density 
ratio of heavy gas contained in the vortex. While the 
specification of w, effectively determines the inner 
boundary of the vortex, the value of r* corresponding 
to the outer boundary is arbitrary. The numerical 
solutions must, however, be terminated for some r* 
greater than unity. In the calculations reported in 
reference 1, they were terminated when the density 
ratio w reached w,,/10 for g, < 2mcerity, and when w 


reached W,,, fOF 2m > Zm(crit)- 


Solutions for w and T* 


Solutions for w and 7*, for two subcritical values of 
the heating parameter, are shown in Figs. 2 and 3. 
The heating rates are 0.17 2n¢mqr) and 0.34 2ncmar). For 
the chosen values of A/,, = 0.7 and wy, = 1.0, Smcerit 
Zm(mar) = 0.49. In this range of 2m / m(maz), the w pro- 
file is rather insensitive to g,,, even though a consider- 
able amount of heating takes place. If g,, were zero, 
T* would decrease quite rapidly with decreasing r* 
because of the expansion in the pressure gradient pro- 
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Fic. 5. Variation of dimensionless temperature, 7*, with dimen- 
sionless radius r*; for g,, greater than g» (crit). 
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Fic. 6. Variation of dimensionless fuel concentration, m.*, with 
dimensionless radius, 7*, for g,, less than g,, (crit 


duced by the vortex. It follows that for sufficiently 
large r*, 7* would increase with r*. In Figs. 2 and 3, 
the solutions are carried out approximately to the 
value of r* at which d7*/dr* = 0. 

As g» increases, d7* dr* increases until the critical 
value, 2micrit), 18 reached. The solutions then change 
character to that shown in Figs. 4 and 5. The essen- 
tial differences are that w has a second stationary point 
for some r* > 1, and that d7*/dr* becomes increas- 
ingly negative with increasing 7*, instead of eventually 
becoming positive, as in Fig. 3. The behavior appears 
to be singular, with w approaching infinity and 7* 
approaching zero as r* approaches some value. Thus, 
in principle, very large temperature ratios, from the 
outside of the vortex to the inside, can be obtained 
with sufficiently large heating rates. 

However, the irregular behavior of w and 7* implies 
very large concentrations of heavy gas in the outer 
portions of the vortex tube. This is illustrated in Figs. 
6, and 7, which give the profiles of m)* corresponding to 
the w and 7™* profiles of Figs. 2 and 3, and 4 and 5, 
respectively. For the subcritical heating rates, the 
profiles of m2* and w are at least qualitatively similar; 
however, for the supercritical heating rates, m.* be- 
comes very large in the outer portions of the tube. 
This is, of course, principally a reflection of the fact 
that the gas density is much larger in the outer portions 
because of the large temperature change. 

Values of 2mccrit) Zm(mar) Obtained from a series of 
integrations are summarized in Table 1. 

The sharpness of the peak in the w distribution in- 
creases rapidly with increasing J/,,,, as shown in Fig. 8 


for three cases having 2»/2m¢crit) ~ 1.02. Thus, as 


TABLE 1 


Critical Heating Rates 





Mim Wa = 0.2 w, = 0.5 ww, =10 w, = 40 
0.5 0.455 0.450 0.441 0.410 
0.7 0.505 0.499 0.487 0.484 
1.0 0.580 0.574 0.557 0.470 
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the Mach number increases, the amount of heavy gas 
which is retained in the vortex, for a fixed value of 
W,, is reduced. On the other hand, increasing .\/;,, in- 
creases the mass flow capacity of the vortex, as may 
be seen from Eq. (5S). 

In contrast to the strong dependence of the profiles 
of won M,,, their dependence on w,, is quite weak 
W/W, is shown in Fig. 9 for three values of w,,, and for 
Lm Zmicrit) = 1.02, as before. This insensitivity of 
WW, to W,, may be understood from a careful examina- 
tion of Eqs. (18) and (19). If we consider only the 
range of w where w (1 + w) > M. (IM, + My), they 


reduce approximately to: 


d log T* 


* 


d log r 


. [me ar w l 
YM im* - 13 — a, oe Ea 
my Cy, (1 + w)? \e* TF 


and 


dlogw — y¥Min” (™ i)( ) 
d log r* T*'2 \m, pe T HL)? 


In the energy equation, the first term represents the 
energy transport by diffusion, and this term depends 
directly on the magnitude of w; however, it is zero at 
two points (the maximum and minimum of w) on the 
curves of Fig. 9, so it seems reasonable that it is fairly 
small over most of the vortex. Of the other two terms, 
the last depends on w w,,, and not on w,, itself. Thus, 
it may be said that the profiles of w w,, are similar for 
different values of w,, because the diffusive energy 
transport is comparatively unimportant 


Performance Parameters 


In addition to the three parameters determining the 
solutions, namely, .\/;,,, w,,, and g,,, the value of r* cor- 
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Fic. 7. Variation of dimensionless fuel concentration, m.*, with 
dimensionless radius, 7*, for g,, greater than g,, (crit) 
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responding to the periphery of the vortex must be 
specified in describing an actual vortex. For given 
values of Mim, Wm, and g,,, the pressure and temperature 
ratios across the vortex, for example, depend on the 
value of r,*. In order to reduce the number of param- 
eters by one, a rather arbitrary choice of r,* will there- 
fore be made. For the cases with subcritical heating 
rates, r,* will be taken as that value of 7* greater than 
unity for which w/w, = 0.1. For the supercritical 
cases r,* will be taken as the value of 7* greater than 
unity for which w = w,,. 

With these choices, there is then, for given values 
of w, and Mi,, a definite relationship between the 
temperature ratio across the vortex, 7./7),, and the 
heating parameter. This relationship is shown in 
Fig. 10 for Mi, = 0.7. It will be noted that the char- 
acteristics are markedly different for the subcritical 
and supercritical solutions. As the heating parameter 
is increased, the temperature ratio, 7\./7,, increases 
continuously until gp cris) is reached. At this point, 
the profile of w switches from the type shown in Fig. 2 
to that shown in Fig. 4, and the temperature ratio 
becomes very large. As g, is increased beyond 
Zerit), 1/1, decreases again, because with the above 
choice for the outer radius, 7,* decreases with increasing 
Lm. FOr 2m/2mimar) = 1, the outer boundary of the 
vortex is at the maximum of w, ie., at r* = 1. Be- 
cause of the sensitivity of the profiles in the neighbor- 
hood of £m crit), it was rather difficult to establish the 
exact end points of the subcritical characteristics; hence, 
these end points are labeled approximate. For MW), = 
0.7, the maximum values of 7°./7;, which can be ob- 
tained with a simple maximum of w vary from about 
3.0 for w,, = 0.2 to 2.5 for w,, = 4.0. Values for other 
Mach numbers are given in Table 2. 

It is convenient to relate other performance param- 
eters to the temperature ratio, w,, and J/,,,. The 
pressure ratio from the inner to the outer boundary 
is shown in Fig. 11 for J, = 0.7. As w,, increases, 
the pressure drop through the vortex increases also. 
This increase may be regarded as due to the increased 
resistance to flow offered by the more dense distribution 
of heavy gas as w,, increases. 
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Fic. 8. Effect of the tangential Mach number at the point of 
maximum density ratio, /,,,, on the variation of the density ratio, 
w, with dimensionless radius, 7*. 
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Fic. 9. Effect of the maximum density ratio, w,, on the var 
iation of the density ratio, w, with the dimensionless radius, 
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The Mach number of the vortex at its periphery, 
M jp, is related to M,,, by 


Mip = Mum/Tp*To*!? (25) 


Since w/w,, and 7* were found to be independent of 
w,, for given .M/,,,, it is not surprising that, as a function 
of 7./T,, Mi, depends only on M,,,._ This dependence 
is shown in Fig. 12, from which it is clear that rather 
high values of /,, are required when 7), 7’, is small. 
This is due simply to the very rapid decrease of 7),* with 
increasing r* for these cases; however, it is important 
to note that the mass flow capacity of the vortex de- 
pends on M,,, [Eq. (S8)], and that as 7, 7. (or T,*) de- 
creases, it will become more difficult to obtain satis- 
factorily high values of J/,,,. 

From Eq. (8), the mass flow capacity of the vortex 
per unit length, 31, is related to \V/,,, and the diffusion 
coefficient. Using Eq. (16), the relation may be written 


as 
E 3[(2rkm,(1 + m,/mz) |'/? [me 
IN, _— rm = l x 
8 dy? my, 
7yI,' 2 
2 Vv (26) 
[1 — (W./Wm) 


where VW = Mi,?/T,*'*. Since 7,* is independent of 
Wm, SO iS V; its dependence on M,,, and 7,7, is shown 
in Fig. 15. It is clear that for given values of 7, and 
Mim, SN increases with increasing 7./7,,; however, if 
M1, is held fixed, Y decreases with increasing 7°, 7, 
as indicated by the lines of constant J/,, included 
in Fig. 13. 

The radius r, of the inner boundary of the vortex 
efiectively determines the maximum possible length-to- 
diameter ratio of a vortex such as that in Fig. 1, since 


TABLE 2 
Maximum of 7./7, for Subcritical Solutions 


M im Wm 0.2-1.0 Wm = 4.0 
0.5 4.2 3.8 
0.7 3.0 2.& 
1.0 » 2.4 
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all of the mass flow must leave the vortex axially near 
its center, and the radius of the exit nozzle must be 
less than r,.. Thus, the ratio r,/r, is an important per- 
formance parameter. Its dependence on M,,,, w,,, and 
7. T, is shown in Fig. 14. In contrast to \/,, and 
VY, 7. rp is dependent to some extent on w,,. This is 
because the profile of w depends on the magnitude of 
w itself, and not w/w,, alone, near 7,, where w ~ w,. 
All the values or r,/7, are near 0.5 or larger, so that for 
the range of M/,,, covered in Fig. 14, the heavy gas is 
essentially confined to the outer half of the vortex 
radius. 

It is convenient to refer the heating parameter to the 
conditions in the vortex at its periphery. Thus, 

») 2 
Z=- Sy = Sm(Mop*fp*?/Tp*) (27) 
M1 + w.)CnT> 

is a direct measure of the rate of enthalpy increase in 
the outer portion of the vortex, and if it is given as a 
function of 7. T,, w,, and A/,,,, it permits the calcula- 
tion of the neutron flux, ¢, required to give any speci- 
fied temperature ratio. The dependence of g, on W,», 
and 7.7, for Mi, = 0.7 is shown in Fig. 15. It in- 
creases rapidly as 7.7, increases, as would be ex- 
pected. The large difference between the supercritical 
and subcritical cases for a given 7’. 7, is due simply 
to the fact that w,/w,, = 1 for the former and w,/w,, = 
0.1 for the latter. If the values of g, for the subcritical 
case are multiplied by 10, they very nearly fall on the 
curve for the supercritical case, in the range of 7. 7, 
where the two overlap. 

Finally, to complete the description of the vortex 
from the nuclear standpoint, we need some measure of 
the total amount of fissionable gas and of propellant 


gas init. We may define a mean propellant density by 
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Fic. 1]. Dependence of the overall pressure ratio, p./p,, on the 
overall temperature ratio, 7/7; for Mj, = 0.7 


hd 
Pice*re*? + 9 | pi *r*dr* 
Jr 


a * ‘i 
- S S) 
*K *? = 
Pip Pip lp 
if we assume p;* is constant at p;,* for r* <r,*. Simi- 
larly, if we define a mean density ratio, ®, by 
oh Tr ¢ *rp* 
u Ps 1 , p ‘ 
=2-? al — { no*r*dr* (29) 
Wy Pi Ppp > oY re* 


then the average density of the fissionable gas is p2* 
@(pi/ Pip) Pip. Pi/ Pip and ®@ w,, are given in Figs. 16 
and 17 for M/,,, = 0.7. We see that ®@/@,, is largest for 
the supercritical cases, because w = w,, both at r* = 1 
and at r* = r,*; and that it has a maximum, because 
w deviates more from w,, in the outer portions of the 
vortex as 7),/7. decreases. &, pip increases mono- 
tonically with 7,/7., simply because of the influence 
of the temperature on the density. There is also some 
influence of w,,, which arises from the variation of the 
pressure drop with w,,. 


Performance and Criticality Estimates 


The significance of the results obtained from the 
preceding analysis of the heating separation process 
may perhaps be better revealed if incorporated into 
the performance analysis of a gaseous vortex reactor 
for rocket propulsion. In the present section the char- 
acteristic performance parameters of the rocket are 
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calculated in terms of the parameters derived from the 
separation analysis and from the nuclear configuration 
To be definitive, such an analysis would have to include 
an optimization of the entire system for some mission, 
and this would imply a detailed study of many aspects 
of the vehicle configuration, as well as those of the 
power plant. The intent of the present analysis, how 
ever, is much more modest; it is simply to provide a 
physical feeling for the several parameters which 
characterize the performance of the vortex tube as ap- 
plied to rocket propulsion. 

To this end, the requirements imposed on the reactor 
system by the overall energy balance, and by criticality, 
will be estimated. The weight, thrust-to-weight ratio, 
neutron flux, and various other parameters will then 
be given for some representative examples. 

It is assumed in these computations that the reactor 
core consists of a cylindrical bundle of vortex tubes 
such as that shown in Fig. 1, the diameter being equal 
to the length / of an individual tube. Both graphite 
and beryllium oxide will be considered as moderators. 
In those cases where the temperature ratio 7. 7), is such 
that heat must be added to the propellant, by fission- 
able material, before it enters the vortex tubes, it will 
be assumed that the fissionable material not contained 
in the tubes is uniformly dispersed in the moderator. 
Some reduction in overall size and weight should be 
obtainable by concentrating this solid fuel in as small 
a volume as possible, but this refinement will not be 
considered here. 


Energy Balance 


Although most of the fission power generated in the 
vortex tubes appears directly as an increase in the 
enthalpy of the gas mixture, there will always be some 
fraction ¢ carried by radiation which penetrates the 
gas and is deposited directly in the surrounding solid 
regions of the reactor. The radiation mechanism is 
both thermal and nuclear (gamma rays and neutrons) ; 
therefore, in general ¢ will depend on the properties of 
the gas and the nuclear characteristics of the reactor. 

In the special case that the propellant gas is the only 
available heat sink, the propellant must be passed first 
through the solid regions of the engine to remove the 
energy deposition represented by ¢, and thence into 
the vortex tubes for further heating. Consider now 
a system in which some of the fissionable material is 
contained in solid fuel plates so that a fraction f/f of 
the fission reactions occur in the solid regions. If P 
is the total power generated in this system, /;, the en- 
thalpy per unit mass of propellant corresponding to the 
temperature in the vortex core, and m the mass flow 
rate of propellant, then P = mh,. The total power 
deposited in the solid regions P, is then 


P,=fP+¢éP.= ([f+¢0 -—f|P 


where P, = (1 — f)P, the fission power released directly 
into the gas mixture in the cavity. Further, if h, is 
the enthalpy corresponding to the maximum tempera- 
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ture in the solid, then for the present case P, = hp, 


and we obtain 
hh, = 1/|f+ ¢1 — f)] (30) 


Now, the specific impulse / for a rocket engine with a 
completely expanded nozzle is simply g*/?/2 = h, where 
g is the gravitational constant. Thus, the ratio given 
in Eq. (30) is the square root of the corresponding spe- 
cific impulses. The significance of this result may 
be appreciated if one thinks of the quantity /, as the 
maximum possible specific impulse attainable in a 
solid-fuel (temperature-limited) reactor. Then /,//, 
would be the performance gain to be realized in in- 
corporating some gas-phase fission heating. Fig. 18 
shows this ratio as a function of f for the case ¢ = 0.1. 
Evidently, with an all gas-phase reactor (f = (), it is 
possible to achieve specific impulses about three times 
as large as with an all solid-fuel nuclear engine. 

If the restriction that all heat must be transferred 
to the propellant is removed by providing radiators for 
rejection of heat from the solid parts of the reactor, 
h. h, can, at least in principle, be larger than the value 
given by Eq. (30). However, if an appreciable frac- 
tion of the total power is radiated, the radiator weight 
will be so large that the thrust-to-weight ratio will 
necessarily be much less than unity. This is not 
necessarily a disadvantage for space propulsion appli- 


cations. 


Criticality 

Estimates of the critical size of the reactor have been 
obtained from two-group, two-region calculations for 
a completely reflected cylinder. The calculations 
were done on a digital computer using a two-group, 
two-region code, and a reflector savings program.* 

The reactor core was taken to be a homogeneous 
mixture of moderator, plutonium and hydrogen. The 
proportions of the latter two were determined from the 
separation analysis. Fast fission was neglected and 
the resonance escape probability was taken as unity. 

The thickness of the beryllium reflector for each 
reactor and the void fraction 8 of the core were selected 
to minimize the combined core and reflector weight. 
The calculations made by Lafyatis and Nelson indi- 
cated that 8 = 0.6 would be near optimum for most 
configurations of interest. 

The fission cross section of plutonium has a marked 
resonance at a neutron energy of 0.3 electron volts, 
which coincides with the mean thermal neutron 
energy (Maxwell distribution) for a temperature of 
2,.330°K. Thus, the neutron temperature has an im- 
portant effect on the critical mass for the systems con- 
sidered here. It was assumed that the neutron tem- 
perature was equal to the moderator temperature, 


which in turn was taken equal to 7;. 
Rocket Engine Performance 
For the present purpose, the performance of the vor- 


* The calculations were made by P. G. Lafyatis and M. L. 
Nelson of the Oak Ridge National Laboratory, USAEC 
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tex tube reactor may be characterized by the specific 
impulse, the total thrust, and the ratio of thrust to 
reactor (core plus reflector) weight. 

The specific impulse is determined from the total 
energy content //, of the propellant gas at the vortex 
“core’’ condition. In the present system, the total 
energy is given by the enthalpy per unit mass of gas, 
h., plus the kinetic energy due to the tangential ve- 
locity v, of the propellant at the core radius 7,.. Thus, 
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H, = h, + (22/2) (31) =. d 
Pt = Pre | d 
If we use the fact that the tangential velocity in the y+ 
vortex is inversely proportional to the radius, then we and , 
can relate v, to vp, the tangential velocity at the tube 
periphery, thus: si =« 3% € =3 h ‘ 
, ; y+ 1 V 
Ue" = V¥eReT pM ip*(1p fo)" t] 
v aren Thus, Eq. (34) may be written FE 
where y,, X, refer to the core conditions; and Eq. (31) 
may be written My! = wry*pic[(y — 1)h-|'? X : 
tf 
alae ‘ y¥+1)/2(y-1) 2 , 
7 tee ] e Vn \~ 7 ) BY; 2 le " Yt ‘ o> E 
H,. = h.| 1+ M,,? | - ~+ (32) ( (35) 
2 Te (ae i ae Ty Fe h 
where we have used h, = C,,.7°... For the systems to be Substitution into the thrust equation yields 
considered here, the second term in Eg. (32) is in the eo A 
order of 0.2 and will be neglected in these computations. F = (w/4)Bl?yp, ( _ ) x 
This means that we are underestimating the energy ale F 
available to the gas and will therefore obtain a slightly ) V+ /20-D fy \2 fp \2 al 
» < c 90 
conservative estimate of engine performance. (. rm +) (: ) (° ) (36) tt 
The total thrust F from the assembly of vortex tubes : 
in the reactor is where we have used the relation p;. = yp./(y — 1)h,- and re 
F = Nigc (33) introduced the expression for the reactor core volume 0 
hile i mr,’ N1/B, which is 7z/* 4 for a square cylinder, to 
where / is the length of the tubes, V the number of tubes, eliminate VN. The quantity 7, 7, is the ratio of nozzle 
IN, is the mass flow rate per unit length of tube, and c throat radius to the radius at which w — w, in the vor- 
is the rocket exhaust velocity, where c?/2 = h,. Now, tex tube. Note that 0 <r, 7. < 1. It can be seen Ww 
the product S1,/ is related to the nozzle throat radius from Eq. (36) that the thrust of this system is actually el 
r, by limited by the size of the area which can be provided 
2 " ~ at the back face of the reactor for the exhaust nozzles, 
M,/ = awri2 pias (34) : ‘i 5 ' 
without allowing the fissionable material to be blown tl} 
where p; is the propellant gas density and a, the ve- out. This area is 7§/*(r,/r,)?/4. The thrust is also 
locity of sound at the throat. These are related to proportional to the pressure p,, but independent of the o 
the core conditions by “core” temperature 7. The ratio 7,7, is of course 
TABLE 3 
Numerical Examples of Vortex Tube Reactors* 
Case Number al 
Parameter ] 2 3 5 6 a 
aK 1.56 1.56 1.56 4.05 4.05 4.05 + 
Wm 0.5 1.0 4.0 0.5 1.0 4.0 7 
Mp 1.0 1.0 1.0 1.0 1.0 1.0 N 
75, ae 7.020 7,020 7,020 10,000 10,000 10,000 
Des 4,500 4,500 4,500 2,470 2.470 2,470 
pe, atm. 100 100 100 100 100 100 h 
Pp, atm. 213 231 498 171 188 338 ; 
Yo) Ko 0.66 0.66 0.65 0.59 0.59 0.58 pr 
W/m 0.34 0.35 0.36 0.56 0.56 0.54 of 
Lp 0.96 0.96 1.50 56 56 81 
IM, Ib. /sec. ft. 0.0198 0.0198 0.0198 0.0121 0.0121 0.0121 m 
I, sec. 1.196 1,196 1,196 1,426 1,426 1,426 Ww 
ng, em. 0.40 XK 1018 0.85 K 1018 5.2 X 10% 0.49 * 10'8 1.01 X* 1038 5.12 X 10% 
Ett: 15.4 10.2 3.8 16.8 11.7 1.5 be 
(14.0) (8.5) (4.1) (16.7) (10.1) (4.5) at 
F/(ri/re)*, Ib. 15.3 X 106 6.7 X 108 0.90 & 108 14.6 X 108 7.02 X 10° 1.01 & 108 ‘ 
(12.6 X 106) (4.65 X 108) (1.02 X 108) (14.3 X 108) (5.25 X 10*) (1.01 X& 108) Cc 
Weight, Ib. 180,000 75,000 10,000 230,000 100,000 15,000 me 
(120,000) (36,000) (5,300) (180,000) (55,000) (6,700) 
[F/(r:/re)®|/weight 85 89 90 64 70 67 ce 
(98) (133) (192) (80) (94) (150) 
N/(rt/tr-)? 55,400 36, 700 13,300 66, 200 46,100 17.100 
(50,400) (30,600) (14,300) (65,800) (39.800) (17,100) tu 
o/(r:/%, 2, ac 
neutrons/cm.?-sec. 0.78 XK 1038 0.53 & 108 0.26 * 1038 1.6 X 1038 1.0 * 10%8 0.54 & 10% 
(0.85 X 1018) (0.64 & 1038) (0.24 X 10'8) (1.6 X& 1038) (1.2 X 10%) (0.54 & 10138) by 
rp(re/Pc), in. 0.305 0.248 0.154 0.304 0.253 0.160 on 
(0.291) (0.226) (0.159) (0.303 ) (0.236) (0.160) : 
VI 
* Performance parameters based on 6 = 0.6, y = 1.31, Q = 7.86 X 107! cal./fission, Cp; = 3.54 cal./gm.-°K., and of = 975 barns cl 
for Tp, = 2,470°R., and Cy, = 5.25 cal./gm.-°K. and o; = 810 barns for T, = 4,500°R. . 
asi —— Sa 
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determined by the diffusion process in the tubes. It 
depends principally upon \/,,,, increasing as M,,, in- 
creases. 

Eq. (36) contains the factor (r;,/7,)?, which may have 
any value from zero to unity. Thus, for a given set of 
vortex tube parameters and critical reactor size, the 
thrust may have any value from zero to that given by 
Eq. (36) for r, 7, = 1. As this ratio increases, the 
number of tubes in the reactor increases, each tube 
decreasing in diameter. This may be seen by equating 
Eqs. (33) and (36) and solving for NV. Then, using 
h, =~ Cy-T., we obtain 


y¥+1 


mBip.y'!* ( 2 ) sid (“ ' (“ r (37) 
od 
49, (RT.)!/2? \y + 1 lp le 


For a given set of vortex tube parameters, NM, and 7, 7p 
are fixed. Then, for fixed p, and 7., the number of 
tubes is proportional to (7; 7,)?. 

The average neutron flux required to give the heat 
release rates implied by these performance parameters 
may be determined from Eq. (27), which gives 


(1 -+ WU, art « My 
d = 2p ) 


2nopo sO ET y? 


where i, = CyT),. The last factor above may be 


eliminated by applying Eq. (35). If we note also that 


Nop = WpPip/M2 = Wy(p-/me2RT,) (Pp Pe) (Te/ Tp) 
then the expression for @ may be written 


gp(1 + w.)meCnl (yRT,)'" 
> = x 
2w,l a ,O 
y+1 


eye) ENE « 
(3S) 
r, 77 % Po] \'p re 


gp increases with increasing 7°..7, (Fig. 15), and C,, 


also increases with increasing 7°, so that ¢ will increase 
quite rapidly with increasing 7’.. 


Numerical Examples 

It was mentioned previously that the enthalpy ratio 
h. hy (or the temperature ratio 7, 7,) is the most im- 
portant parameter in determining the characteristics 
of the vortex tube. On the other hand, the fissionable 
material concentration determines the overall size and 
weight of the reactor. Accordingly, examples have 
been selected with two representative values of 7. 7,, 
and three values of w,,. The pressure at the vortex 
core has been taken as 100 atmospheres for all cases, 
so that w,, is nearly a direct measure of the fuel con- 
centration in the gas phase. 

Cases 1, 2, and 3 have 7. 7, = 1.56. This tempera- 
ture ratio corresponds to about half the heat being 
added to the propellant by solid fuel elements, and half 
by the vortex tubes. The temperature leaving the 
solid reactor, 7,, has been set at the upper limit for 
graphite fuel elements, namely, 4,500°R.; hence the 
chamber temperature, 7,, is 7,020°R. This is of the 
same order as the chamber temperature for contem- 
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Fic. 17. Dependence of the heating parameter, g,, on the overall 
temperature ratio, 7,,/7 for Mim = 0.7 


porary chemical rocket motors—e.g., liquid O, and 
JP-4 give T, = 6,300°R. 

Cases 4, 5, and 6 have the largest values of 7, 7, pos- 
sible for 7, = 10,000°R. They represent systems with 
chamber temperatures considerably higher than the 
best obtainable with chemical rockets. At the same 
time, the temperatures entering the vortex tubes are 
considerably below the limits for graphite or the refrac- 
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Fic. 18. Specific impulse ratio, J-/Ip, as function of gas-phase 
fission heating. 


tory metals. These systems have all fissionable ma- 
terial in the gas phase. 

Some characteristics of these reactors are shown in 
Table 3. The values for beryllium-oxide-moderated 
reactors are shown in p :entheses. The first few 
quantities listed are charac eristic of the vortex tube 
itself, and hence are independent of the moderator 
used. 

Although the reactor weights (core plus reflector) 
are very high for fuel concentrations of the order of 
0.5 X 10% cm.~*, they are quite reasonable if the con- 
centration can be increased to about 5.0 & 10'8 em.~* 
These higher fuel concentrations imply high pumping 
pressures, however. For Case 3, p, is 498 atmospheres, 
or 7,320 psi. The pumping pressure must be taken as 
about twice this figure, to allow for the pressure drops 
in the inlet nozzles and solid reactor. It should be 
noted that the increase in average fuel concentration, 
nz, by a factor of 10 from Case 1 to Case 3 requires only 
a little more than doubling of the pumping pressure. 
This increase is due to the higher pressure drop caused 
by increasing w,,. To achieve the same concentration 
increase by raising the general pressure level with con- 
stant w,,, it would be necessary to increase the pumping 
pressure tenfold, to about 4,000 atmospheres. The 
beryllium-oxide-moderated reactors are in all cases 
lighter than the corresponding graphite-moderated 
reactors. This is due to the decreased neutron-slowing- 
down length and decreased diffusion length in beryllium 
oxide as compared to graphite. 

Although the reactor weights are rather high, the 
thrust-to-weight ratios are also quite high if r,/r, is 
near unity. It may, however, be more realistic to 
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take r, r, a little less than unity, to allow for three 
dimensional flow effects in the long, thin, vortex tubes 
If r,/r, is taken as 0.5, for example, the beryllium-oxide 
moderated reactor described in Case 3 has a thrust 
of 255,000 Ibs., and a thrust-to-weight ratio of 48. The 
latter figure is somewhat lower than that for a chemical 
system; however, the specific impulse is at least twice 
that of the best chemical rockets. In order to give 
this performance the reactor would contain 950 vortex 
tubes, each 0.32 in. in diameter. The average neutron 
flux would be 5.2 X 10'® neutrons ‘sec.-cm.? 

It seems from comparison of Cases 3 and 6, for ex- 
ample, that as 7, is increased, the critical size and 
weight of the reactor increase. This effect is due to 
a lower average fuel concentration in the higher tem- 
perature reactor. The gas-phase fuel concentrations, 
f»,, are about the same in the two cases, but the lower- 


“9 
temperature reactor contains considerable solid fuel. 


Concluding Remarks 


In conclusion, it seems worth while to summarize 
the more important results of the analysis. One of 
these is that the mass flow capacity per unit length of 
the vortex tube, as limited by the rate of diffusion of 
propellant through the fissionable gas, is independent 
of the tube diameter. Thus, even though the mass flow 
per unit of tube length is only of the order of 0.01 pounds 
per second per foot, the mass flow capacity per unit 
volume of an assembly of small-diameter tubes can be 
quite large. 

A second important conclusion is that containment 
of the fissionable material within the vortex, away from 
its outer boundary, is possible only if the pressure 
gradient decreases monotonically with increasing ra- 
dius. Asa corollary of this, there is a critical heating 
rate above which such containment is impossible. 

Insofar as the vortex containment mechanism limits 
the performance of the gaseous nuclear rocket, the 
performance capability is very high. It may well be 
that the performance will actually be limited by other 
factors, such as radiant heat transfer from the gas to 
solids, the difficulty of maintaining fissionable material 
in gaseous form at the required concentrations and 
temperatures, and the difficulty of generating vortices 
with the required low radial mass-flow rates and high 
tangential velocities. The potential of the gaseous fis- 
sion rocket seems to warrant intensive study of these 


problems. 
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Torsional Analysis of Thin-Walled Shells 
With Many Closely Spaced Spar Webs 


DAVID W. FOX* anp ROBERT M. RIVELLO** 
The Johns Hopkins Unwersity 


Summary 


The torsional shear flow distribution and rigidity of thin-walled 
cylinders with closely spaced internal spar webs is investigated. 
The spar webs are assumed to be so numerous that the spar 
as a continuous shear-resisting me- 
conditions 


material may be regarded 
dium. The differential 

which must be satisfied by the shear flow are derived from the 
This differential equation is solved 


equation and boundary 
Bredt theory of torsion 
explicitly for those cases in which the shell cross section is given 
by a parallelogram, trapezoid, triangle, or double wedge. The 
results are related to the corresponding rigidities and shear 
flows given for the same section by the Bredt theory in which 


the internal spar material is neglected. 


Symbols 
A area 
G = shear modulus of elasticity 
J = torsion constant 
1 torque 
a,b,c £ . : 
i. 4" cross-sectional dimensions 
h = height 
k = core solidity defined by Eq. (6) 
n = number of cells 
q = shear flow 
s = arc length 
t = thickness 
x,y = Cartesian coordinates 
a = angle between section face and x axis 
0 = ‘angle of twist per unit length 
T = shear stress 
Subscripts 
i refers to ith cell 
t top 
b bottom 


Introduction 


— CURRENT TREND toward increased unit loading 
and decreased thickness of aerodynamic surfaces 
has resulted in the use of structures having relatively 
heavy skins supported by a large number of internal 
webs that stabilize the skin and resist shear and tor- 
sional loads. For this type of structure, known as 
multispar or multicell construction, the torsion may be 
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analyzed, when warping is not restrained, by the well- 
known Bredt theory.!' When both the shear flow and 
the angle of twist are desired this procedure requires 
the solution of a system of simultaneous linear alge- 
braic equations equal in size to the number of cells plus 
one. When the number of cells is large and exact 
solution of the simultaneous equations becomes tedious, 
approximate methods of solution involving less labor 
become desirable. In this direction, methods of suc- 
cessive approximation and iteration?‘ have been used. 
The process of idealizing the structure by lumping 
several webs together in order to reduce the number of 
equations to a more manageable size seems less desir- 
able. On the other hand, if the number of webs is 
quite large, another approach suggests itself, the con- 
cept of a continuous shear-supporting spar web medium. 
In this idealization, which is the subject of this paper, 
the simultaneous equations are approximated by a 
linear second-order differential equation with end 
conditions. In this note the boundary value problems 
are solved explicitly for the differential equations that 
arise from the consideration of cross sections in the 
shape of a parallelogram, trapezoid, triangle, or double 
wedge. These results, which apply to sections often 
used in supersonic airfoils, should be useful for design 
studies in which the work involved in solving a large 
system of simultaneous equations is not justified or when 
the number of cells is so large that solution of the sys- 
tem of equations is not feasible. An example of the 
latter case occurs in the consideration of some forms of 
sandwich box-beam construction.® 

The theory presented here was given in less complete 
form in reference 6 and has been extended to the cal- 
culation of shear flows due to a shearing force and to 
the location of the shear center.’ 

A recent note by Niedenfuhr* treats the problem of 
torsion of a uniform rectangular section with any num- 
ber of equally spaced identical cells by the techniques 
of linear difference equations with constant coefficients. 

In contrast to this paper, in which it is assumed that 
the continuous shear-supporting web medium possesses 
shear rigidity in one direction only, Seide® and Cheng" 
have treated the case of a rectangular cross-section 
sandwich with a nonisotropic core and McComb!! 
studied rectangular, triangular, and diamond cross 
sections with an isotropic core reinforcement. 


Theory 


We consider a cylindrical structure with the cross 
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Fic. 1. Segment of multicell structure 


section divided into a number of cells by vertical in- 
ternal webs. We assume that all walls of the cross 
section are thin and that warping is unrestrained so 
that Bredt’s theory of torsion applies. 

For such a structure the 7th cell together with the 
two adjacent cells is shown in Fig. 1. In terms of the 
shear flows, the angle of twist per unit length is deter- 
mined by the Bredt equation! which in the usual nota- 
tion is given by 


2A 0 = [. (q/Gt) ds (1) 


Writing out the integral and the expression for .1; in 


Eq. (1) gives 


saeaad dy 
(Ji — 91-1) F — (Ji41 = ne) x 


nb (xi to(y) 
‘e dy 4 [- jl E 1 (2)") + 
J hn(xi) ty(y) 4 J Xi-1 Lé(x) dx 
] dh, 2 ade a ‘ 
1+ > dx = 2G0 h(x)dx (2) 
t,(x) dx f a 
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in which h is defined by A(x) = h(x) — h(x), the sub- 
scripts ¢ and 6 referring to the top and bottom surfaces, 
respectively. Here we have assumed G to be a con- 
stant throughout the structure. If G is not the same 
everywhere, Eq. (2) still holds if the G occurring there 
and in later equations is interpreted as a reference 
shear modulus G,,; and the ?’s are interpreted to be 
equivalent thicknesses defined by 


t = t(G/Grer) 


equiv 


Assuming that h,, fy, t,, and t, are smooth enough, the 
mean value theorem for integrals may be applied to 
Eq. (2) to yield 


) h(x;-1) h(x;) 
te “ Ge. 5s ote J i4 - = 
, ~ b(X4-1) ita tie Ealx,) 
B(&)gi(xi — Xi-1) = 2GOh(z) (x; — Xin) (3) 


where Xsatisfies x,) < * < xj, f,(x,) is the effective 


thickness of the cell wall at x; defined by the relation 


l I ms) dy 
. = - (4 
bey(X3) h(x;) Jh(xi) ty(y) 


and @(x) is given by 
Bx) = 1 E ‘ (*)] 2 ‘ 1 E 1 Gall 2 
oti t,(x) dx th(x) an7 3 


Define k(x;) by the relation 


R(x) = bel(xi)/(Xin1 — Xi) = be(1)/bw; (6 


By use of this definition Eq. (3) can be put into the sug- 


gestive form 


(Gita — Gi) h(x;) (gi — Gi-s) A(xi-1) 


(Xin. — X;) R(x;) (x; — Xi-1) R(xi-1) 
Xi X41 
B(a)g; = —2G0h(z) (7 


which is a formal difference approximation to the 
second-order linear differential equation 


d (7 dg 


) — B(x)g = —2G6h(x) (8) 
dx \k(x) dx 


The boundary conditions that are appropriate with 
Eq. (8) arise from the consideration of the Bredt 
equations for the end cells, and consequently they de- 
pend on the manner in which the section is closed. 
Three typical cases are shown in Fig. 2. 

If the section is closed by a sharp edge as shown in 
Case I of Fig. 2, the appropriate equation for the first 
cell, obtained in the same manner as Eq. (7), is given by 


(2 —_ a) (Re = nt) 
Seu 2 X1 — Xo 
R(x) B(#)g1 = —2GOh(z)R(x1) (9) 


This may be regarded as a formal difference approxi- 
mation to the differential boundary condition 


dg dh ees 
—fk = 6 (10 
E dx 20 = 
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THIN-WALLED SHELLS 


The notation a+ means the limit from the right of a, 
and similarly a~ means the limit from the left of a. 

If the section is closed by a rounded edge as in Case 
II, then the first cell equation may be written 


ie — «€ tia = 
(‘: “) _ as(a)e(e)( 2 ) - 
Xo — xy h(x) 
ee x1 — Xo 
= 260 x) ( ) (11) 
h(x) 


and the corresponding boundary condition for Eq. (8) 
is given by 


FE e(! h, 4 | P| 
dx 9 t. h hs Ses 
In Case III the closure is accomplished by a web of 


thickness f,. 
be written 


g—-aa\,, . “h(xa) 
h(x.) — R(x1) 91 = + B(#) (x1 — Xe) s 
ie = Dh t(Xo) 


—2G6h(#)k(x1)(x%1 — xo) (13) 


(12) 


Here the equation for the first cell may 


and the appropriate boundary condition for Eq. (8) is 


ig k 
| 3 —- ‘| “% (14) 
dx -: 


Similarly, the boundary conditions for sharp-edge, 
rounded-edge, and web closures at the right-hand side 
of the cross section are found to be 


dg dh ae 
- + kBq = 0 (15) 
dx dx a 
FE + e(? h, + 1 r)| , (16) 
= j 
dx g t,h Z@7 1... 
and | + k a| = () (17) 
dx a 


Whenever the section is interrupted by a thick web or 
a change in construction which reflects itself in a dis- 
continuity in the coefficients of Eq. (8), it is necessary 
to have two conditions in order to join the values of g 
which hold to the right and to the left. If there is a 
web of thickness ¢, at x = d, then the appropriate condi- 
tions arise from the limiting forms of the equations 
which hold for the cells adjacent to the web. They are 


ld 1+) — g(d 
dq _ ge g(d~) (18) 
k dx'x=a+ ta 
ld i+) — g(d-) , 
and “4 = al AC (19) 
k dx',=a- ta 
If there occurs at x = e a discontinuity in construc- 


tion with no thick web, then the shear flow g and (h/k) 
(dg/dx) must be continuous. The matching conditions 


are thus 


h dq _ hdq 


= (20) 
adStivne* 2. 
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Fic. 3. Parallelogram cross-section. 
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Fic. +. Shear flow distribution for typical parallelogram cross 
section with identical end webs 


In the case in which the section is symmetrical about 
a vertical centerline (supposedly at x = s), we need 
consider only half of the section with the symmetric 
boundary condition 
de 
y = (21) 
dx , 
The shear stress +r applied by the faces to the spar 
web medium is given by 


rt = dg/dx (22) 


and an approximation to the web shear flow in the 
finite case is obtained by multiplying this stress by the 
web spacing 6,,. 


1.0 





0.9 
0.8 
= 
2 
= 0.7 
i) 
> 
0.6 
0.5 
0.4 
ie) 10 20 30 40 50 60 70 
rb/2 
Fic. 5. Jsneii/J for parallelogram cross section with identical 


end webs 















































728 JOURNAL OF THE AEROSPACE 

1.0 . ‘ 
0 100 

3 
09 

£ a 

> 2 

= bt,X 

e bo 

$08 |—"“2: — 

7 | 
0.7 | 


20 30 40 50 60 70 
db/2 

FIG. 6. Qsheit/Qmaz for parallelogram cross section with identical 

end webs 


° 
ro) 











wl) NON 
NX 






























































$ 
0.2 
= + 
0.1 
btgr s 
| 
O it 
10) 10 20 30 40 50 60 70 
db/2 
FIG. 7. Qnin/Qsneit for parallelogram cross section with identical 


end webs. 


The differential equation (8) and an appropriate pair 
of boundary conditions form a self-adjoint system for 
which the existence of a unique solution can be guaran- 
teed under reasonable restrictions on the coefficients 
of the differential equation (8). Such a solution gives 
an approximate solution to the interior cell equations 
and the end cell equations which arise from the Bredt 
theory, provided that k(x) is chosen to be a smooth 
function that takes the values k(x;) defined by Eq. 
(6) at the points x; The function k(x) accounts for 
the variation in spacing and effective thickness of the 
vertical webs and is, in fact, an expression for the solid- 
ity of the core of the idealized structure. 

Altogether, the differential equation (8) with its 
boundary conditions defines an idealized model of the 
given structure in which the shape defined by h, and 
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h, and the face thicknesses ¢, and ¢, is preserved, but 
in which the discrete internal webs are replaced by a 
continuous shear-supporting core. 

It is clear that the solution of the Bredt equations 
(1) and that of the corresponding boundary value prob- 
lem will depend linearly on the prescribed angle of twist 
per unit length. Since the torque rather than the 
angle of twist is usually given, the torque may be re 
lated to the angle of twist by the second Bredt rela 
tionship! 


fu stE As (23) 
i=1 


where the summation extends over all cells. For our 
idealized sttucture, the limiting form 


T = 2S gh dx (24) 


is used; here the integral is taken over the entire width 
of the cross section. In either case the torsion constant 
J is defined by the relationship 


J = 1 /G6 (25) 

It is not easy as a rule to obtain explicit solutions of 
Eq. (8) with appropriate boundary conditions; how- 
ever, simple solutions can be found when h, and /, are 
linear in x, & is constant, and where ¢, and ¢, are con- 
stant but not necessarily equal. The simplest and most 
important cases in which k is constant arise from equally 
spaced webs of constant thickness. Within these re- 
strictions, solutions for structures having cross sections 
in the shapes of parallelograms, trapezoids, triangles, 
and double wedges are given in the following. The 
usefulness of these solutions is increased by the fre- 
quent occurrence of such cross sections in supersonic 
airfoils. 

For thin airfoils with sharp or web-closed edges 
described by the boundary conditions for Cases I or 
III of Fig. 2, the terms (dh,/dx)* and (dh,/dx)? will be 
small compared to unity. If, in addition, the top and 
bottom skin thicknesses are constant, 8 becomes nearly 
constant and approximately equal to (¢, + ¢,) / tip. 


Applications 


In this section we will assume that the quantities 
t,, t}, and k are constants and that /, and h, are, at most, 
linear functions of x. 


(1) Parallelogram Cross Section 


When the cross section is a parallelogram, h is con- 
stant and h, and /, are linear with slope tan a (see Fig. 
3). Accordingly, Eq. (8S) may be put into the es- 
pecially simple form 

(d°g/dx?) — \*°q = —2G6k (26) 


in which X is given by 


1/2 
i | oe + 2) (27) 
ht t, cos a 
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THIN-WALLED SHE 
The appropriate boundary conditions for end closures 
of thickness ¢, and /, at the ends x = +6/2 are given by 
Eqs. (14) and (17) so that 


ie a k ‘| 0 (28) 

dx {. x b 2 

2 4 k | - 
dx i. q x +b/2 


The solution of the system, given by Eqs. (26), (28), 


and (29), is 


and (29) 


g = 2kGOX~*(1 — Csinh Ax — Decosh Ax) (30) 


in which C and D are expressed by 


; A(t, — ¢t.) sinh (Ab, 2) 
RL + (A%ate/R?)] sinh Xb + A(t, + t-) cosh Ab 
(31) 
and 
D = A(t, + t.) cosh (Ab/2) + 2k sinh (Ab/2) 


k?)] sinh Ab + A(t, + ¢,) cosh Ab 


(32) 


R{1 + (At, t, 


Application of formulas (24) and (25) to the right- 


LLS WITH 
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and the expression for the shear flow, given in terms of 


the torque, yields 
g = 27k XA? J 


1(1 — Csinh Ax — D cosh Ax) (34) 


In the special case where /, =/,, the expressions are 


Corresponding to Eqs. 


sinh (Ab >) 


(35) 


somewhat simpler. (30), (33), 


and (34), we have 


g = 2kGOr~? (; ~ 
cosh (Ab 2 


cosh Ax 
>) + (At, Rk) 


2) + (bt,r 


aA7/2k) — ‘) 
2) + (dt,A? 2k) 


36) 


J= 


(Ab 2) coth (Ab 
tkhby\~? 
(Ab /2) coth (Ab 


and 


(Ab/2) coth (Ab/2) + (bt,rA7/ 2k) — 
 é (A6/2)(sinh Xb6/2)~—! cosh Ax 


ie 2bh \(Ab/2) coth (Ab/2) + (bt,A?/ 2k) — 1 


It is interesting to compare the results in this simple 
case with those obtained from the Bredt theory when 
the internal webs are completely neglected. In the 
latter case, which we shall designate as the shell case, 
the results are 





















































hand side of Eq. (30) yields the torsion constant in the Jone = +1 -/2bh (38) 
— 4(bh)? 
: and Juan = (39) 
J = shkX~*[b — 2Dd— sinh (AB /2) } (33) (b/cos a) [(1/t,) + (1/t)] + (2h/t,) 
1.10 T | ] 
Pe Er ie eee ae 
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Fic. 8. Comparison of finite number of cells to theory for a typical parallelogram, cross section. 
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so that the ratios of corresponding quantities are given by se 
ay : di 
gq (Ab/2) coth (Ab/2) + (btar?/2k) — (Ab/2) [cosh Ax/sinh (Xb /2) | (a 
== } el 
shell (Ab/2) coth (Ab/2) + (bt,A?/2k) — 1 fio 
and ne 
ST 
J ( in 2k ) (Ab/2) coth (Ab/2) + (bt,A?/2k) — 1 (4) 
J shell bt,r? (Ab/2) coth (Ab/2) + (bt,r?/2k) (3) 
Examination of the expression for g shows that the A comparison of the solutions for sections having a a 
maximum shear flow occurs at the center and that the finite number of cells with the continuous solution ob- ax 
minimum ocurs at the ends. These values are given by tained from Eqs. (41), (42), and (43) is shown in 
' Fig. 8. In the results given there for typical section 
Gmaz \b/2) tanh (Ab/4) (bt,A\?/2k) ; i ‘ 
Jmaz _ \ 3 : (42) parameters, the web thickness /, is decreased as the ed) 
9 . y) 2/9 es p a x 
Gshell (Ab/2) coth (Nb/2) + (btad*/ 2k) | number of cells is increased so that k remains constant. 
and It is evident that as the number of cells increases the 
ratiOS Qmar/Yshett ANd J/JSshe11 approach those of the 
Ymin (bt,r” 2k) . { ~ 6 ‘ ani r thile a a0 F 
= (43) continuous solution rapidly, while the approach of 
Yshell (Ab/2) coth (Ab/2) + (bt,rd?/2k) — 1 Jmin/Ysvert iS Slower. The figure also indicates the 
A graph of a typical shear flow distribution obtained degree of approximation that may be expected when | 
from Eq. (40) en given in Fig. 4. Similarly, Figs. 5 sections with a large number of cells are approximated 
6, and 7 give graphs of the functions in Eqs. (41), (42). by lumping webs together instead of replacing them 
and (43) for a range of values of the parameters which bya continuous spar web medium. ee | 
encompasses most practical cases. It may be observed The shear flow distribution sor @ 6-cell ene | 
that except for very small values of bf,\2/2k which do compared with the distribution obtained from Eq. 
at e) y sm: alues g*/2 sg 
not ordinarily occur in practice, the values of the maxi- (40) in Fig. 4. 
mau shear flow and torsion sepa are affected (2) Triangular Cross Section 
little by the interior structure. For the minimum é ee 
aide Cees iin in aaah tht tie In the triangular section it is convenient to choose 11. 
: the origin as vertex (see Fig. 9) and to designate the col 
height of the section at the thick end by /,. The height ust 
of the section is given by 








for 


y cs | 
— T | h = h.(x/b) (44) the 


| 
| 
| 7 | The top and bottom faces make angles a, and a, with 
| . e ° = 
———)| | the x axis. For this cross section, Eq. (8) may be put 
| . P = 
= —acn Te | in the form Th 

























































































- | (d/dx) [x(dq/dx)| — 67 = —2kG6x (45) sata 
i . oe J 
a ee ee ee ax eal in which 6 is given by ; 
L — Session iz 5 Rb t, cos a, + ty COS ay, (46 
‘ pee ; = ) 
Fic. 9. Triangular cross section. he  tyly COS a; COS a 
The applicable boundary conditions are given by Eqs. cal 
1.4 — 7 (9) and (17), respectively (by putting a = 0, c = 5b). an 
For these, Eq. (45) has the solution ol 
1.2 — 1 ; 
g = GO[((2k/6)(1 + 6x) — Chy(2V dx) | (47) | 
1.0 ee mailed dot 
a - In this equation and those which follow /y and Kv are sal 
2 08 4 SS modified Bessel functions of the first and second kind, Th 
£ 
ua /) J/Jsheliz.20 | respectively, and of order vy (see reference 12). is | 
g : 
7 06 b= 5.38 in. te 
— h= 0.987 in. C (2k/6)(1 + kb/t. + k/t,6) of 
k=4.25x10° = - ee" (48) 0 
0.4 <0] 3 o50)' (k/te)1y(2V 6b) + V8/ b1,(2 Vb) 
a ee §=3.61 in. 
. — aia , (4) 
0.2 Application of Eqs. (24) and (25) to the expression 
fo 43 given by Eq. (47) yields the following expression for the 
°o ol 02 03 04 05 06 O07 08 09 10 torsion constant : mie 
b orl 
. . ee J = (2he/8){ (kb/8)(1 + 2586/3) — aa 
Fic. 10. Shear flow distribution for typical triangular cross — pee 
section. C[V 8b 11(2V'86) — 1.(2V'86)]} (49) 
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A typical shear flow distribution for a triangular 
section is shown in Fig. 10; it has been normalized by 
dividing it by the shear flow obtained when the in- 
terior material is neglected. As may be seen from this 
figure substantial errors in shear flow result from 
neglecting the interior. The core increases the tor- 
sional rigidity by 20 percent. 

(3) Double-Wedge Cross Section 

The differential equation (45) applies to each side of 
a double-wedge section symmetrical about a vertical 
axis provided the coordinates are those shown in Fig. 








Fic. 11 Double-wedge cross section 


11. Since the section is symmetrical, the boundary 
condition at x = 6 is that given by Eq. (21), and, as 
usual, Eq. (10) holds at x = 0. 

The solution for each half of the wedge has the same 
form as that for the triangle—i.e., Eq. (47)—but here 
the constant C has the value 


C = 2kV 6b /[6214(2V 6b) | (50) 


The torsional constant, obtained in the usual way, is 


found to be 


J = (4h.’5)} (kb/6)(1 + 26/3) — 
2k V 6b /[62,(2V 6b) | X 
IV 86 1,(2V 6b) — Io(2V8b)]} (51) 


The solution given here is not limited strictly to physi- 
cally symmetric double wedges, but in fact applies to 
any double wedge which can be put into mathemati- 
cally symmetrical form. 

Fig. 12 shows the distribution of ¢/gshne1i for a typical 
double wedge made up of two triangles, each with the 
same section properties as that to which Fig. 10 applies. 
The result applies whether or not a thick central web 
is present. It is to be noted that here the torsion con- 
stant is increased by 27 percent due to the presence 


of the interior web material. 


(4) Trapezoidal Cross Section 

For trapezoidal cross sections the differential equa- 
tion has the form given in Eq. (45) provided that the 
origin is chosen to be at the projected vertex of the 
section (Fig. 13). The appropriate boundary condi- 


tions are expressed by Eqs. (14) and (17). From 
these g is obtained in the form 


g = 2kGO|(1/82)(1 + dx) — Chy(2V dx) — 
DK,(2V éx)| (52) 


where C and D are given by 


u(a)r(c) — ulc)r(a u(c)p(a) — ula)p(c) 


= ) 


Pla)r(c) — ple)r(a) pPla)r(c) — r(a)plc) 


(53) 
in which we have put 


p(x) = y(x)Io(2V dx) + V 6/x1,(2V bx 
r(x) = y(x)Ko(2 V 6x) — Vb/x Ky(2 V bx) (54) 
u(x) = (1/67) [6 + y(x) + bxy(x) | | 


ind (x) f—kita, #8 ot 55 
an y(x) =< ‘ (5: 
‘ i+ /t., x = cf ak 


The torsional rigidity found from Eqs. (24) and (25) is 


v(c) — va) | (56) 


J = (2kh, cé") 





where 
v(x) = x* + (26/3)x* — 
26x} CIV dx1\(2V bx) — 12(2V bx)] — 
D[V 6xK\(2V 6x) + K2(2V dx) }{ (57) 
The shear flow is given in terms of the torque 7° by 


g = (2k7/J)((1/8)(1 + 6x) — Ch(2V dx) — 
DK(2V 6x)] (58) 


(Continued on page 752) 
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Fic. 12. Shear flow distribution for typical double wedge 
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Tic. 13. Trapezoidal cross section 








Frequency-Response Functions of 
Orthotropic Sandwich Plates’ 


M. P. BIENIEK* ann A. M. FREUDENTHAL* 


Columbia University 


Summary 


A method is presented for the determination of the frequency- 
response functions of the components of deformation and of 
stress in orthotropic sandwich plates. It applies to the case of 
simply supported rectangular plates, loaded by pressure normal 
to their planes. The analysis takes into account the transverse 
shear deformation of the core and the material damping of the 
core and of the facings. The results are presented in the form of 
expressions suitable for numerical evaluation. 


lustrates the application of this theory. 


An example il- 


Introduction 


iy ORDER TO PREDICT the behavior of a structural ele- 
ment subject to the action of random loading (for 
instance, random noise), it is necessary to determine 
the frequency-response functions of various quantities 
such as the components of deformation, the compo- 
nents of stress,etc. This paperdeals with the frequency- 
response functions for a rectangular sandwich plate. 
The loading is assumed in the form of a normal pres- 
sure, uniformly distributed over the surface of the plate 
and random in time. The frequency-response function 
of any quantity S (deformation, stress) is the function 
S(iw), such that the quantity S corresponding to the 
loading g = le’ is represented in the form 


S = S(t) = S(iw)e™ (1) 


The dimensions of the plate and the coordinate system 
are shown in Fig. 1. The analysis takes into account 
the effect of transverse shear deformation of the core. 
The thickness of the facings is assumed to be small in 
comparison to the thickness of the core. 

The basic theory of such plates has been developed 
in two papers by R. D. Mindlin! * in which the detailed 
discussion of the assumptions is presented. 

The materials of the core and of the facings are as- 
sumed to be orthotropic. The elastic constants of the 
core E,,, Ey,, Ez,, Gz, G, G; are defined by the stress- 
strain relations 


Oye = WZGr€yz, Ore = WGyerzy Ory = 2WGrexy l (2) 
Cee = Estee + Esytyy, yy = Egytyy + Exytec) 


and the elastic constants of the facings E’,,, E’y,, Ezy, 
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G’, by the relations 


, y , ‘7 / ) 
Fxg Ee ek ee I en as t 


, or 
CF xy = 2G sb, 
«= Bata + ios 


9 
(.3) 


The transverse shear deformations of the facings are 
neglected. If the materials of the core and of the 


facings are isotropic, then 
G, = G, = G. = G = E/2(1 + ») 
E., = E, = E/(1 — »*), Ey = Ev/(1 — v*) 
with similar equalities for the constants of the facings. 


The effect of damping is introduced by the use of 
the complex moduli of elasticity (reference 3, p. 276) 


I 


G, = G,(1 + inz), G, = G,(1 + in,) 
G, = G1 + in,), Ee. = Ex(1 + ins) (4) 
Ex(1 + ty) 


Ey = E,,(1 + ty), E, = 
The coefficients yu,, uy, wu, are 27 times the specific 
damping AW./W,, where AW, is the energy dissipated 
during one cycle in shear, and W, is maximum shear 
strainenergy. The coefficients 7 are related in a similar 

ray to the dissipated and the strain energies in the 
cyclic tension-compression tests. In general, u;,.. .,n2 
depend on the frequencies and amplitudes of cyclic 
stresses or deformations. In first approximation, it is 
assumed that the damping is frequency-independent. 
This assumption is in fair agreement with experimental 
results for various materials. It is also assumed that 
damping is independent of amplitude of applied stress. 
This assumption is less justified by the experimental 
evidence, since the damping of real materials usually 
increases with increasing amplitudes. However, with- 
out this assumption nonlinear equations of the problem 
would be obtained. A proper selection of constant 
coefficients uz,..., Nz, for an expected (not too wide) 
stress and strain range, may provide sufficient accuracy 
of the results of the analysis. The coefficients of 
damping and the complex moduli of elasticity of the 
facings will be denoted by y’,, ’::, yy n’r, and G’.. 
E’ yxy Ey, E's, respectively. 


Basic Equations 


The equations of motion of an elastic, orthotropic, 
sandwich plate are derived in reference 2, in the form 
OM, , OM, er Oy, 

a 5 
Ox oy ot? 


















































kb s -| where 
x f facing K, = «°G_h, K, = «*G,h 
=— "y 5 ») Pad 2 ) 
: an Dp, = Eh 12 + E’ zzh’f 
D, = E,,h*/12 + £',,h?f/2 (S) 
facing Dy = E,h*/12 + E',,h*f 2| 
f H,, = G,h?/12 + G'h*f/2 
y z 
Fic. 1. Coordinate system and dimensions of plate 1 ~e0,000 
ma 
- $00) 
(3) = \w " y- w/ wD, | 
xy, 4X an 
» are | 
the l> 2 
300} 
Fic. 2. Components of deformation w, yz, wy. 
| 
’ 200 
oM,, Ol, 0 I O*y, | 
Ox Ov . or” 
‘ 100} 
ngs. O00, 4 o?, 4 vu Ow | 
e ol gz j = (5 | 
: Ox Ov d ot? 0} Mima nas ——— wii, 
6) T Ol 02 03 
where M,, M,, Al, are bending and torsional moments, Hen & Deities, dj os tein ol eo 
respectively; @Q,, Q, are transverse shear forces; w, 
(4) ¥,, Wy are the components of deformation of the plate = 
as shown in Fig. 2; g is the loading perpendicular to the a 
plane of the plate; / is the moment of inertia of the 
“ific mass of the plate with respect to its middle plane; 1/ | 
ted is the mass of the plate. The internal forces .\/,,..., | | 
ear Q, are taken per unit length of the cross section, the TT 
° ° ° ° i 
ilar loading g, moment / and mass / per unit area of the * qe SS _oe 
the plate. 
fies If the densities of the core and the facings are p and 
clic p’, respectively, then 
t is — 
nt. M = ph + =e} pee (6) 
tal IT = ph*/12 + p’fh?/2) 
nat ai a a a a se the thickness 
ibe where / is the thickness of the core and / is the thickness Fro. 4. Deflection &, (we = 1, © = 1) 28.0 function of frequency 
“asi of one facing. w/w (plate without damping) 
~t The relations between the internal forces V/,,..., 
y ; : asf waa — pm - . mikes Say Sas 
+h Q, and the components of deformation w, ¥,;, ,, as The equations (5), (7), and (8) are also valid for 
1- _ . c . ° am © , . . . “fs . . . 
derived in reference 2, are plates with damping if the moduli £,,,...,G, are 
2m a aed 
: replaced by the complex moduli £,,,...,G., E’,, 
int QO, = K,[¥, + (Ow/dx)] | =, Skene P : F ~~ 
ie) 0 K, |v, + (dw/dy)] ..., G’,. As it follows from reference 2, the coefficient 
: » = Kyl, w/ Oy , : 
= « for sandwich plates is very close to 1, and therefore 
cy M, = D,(Op,/0x) + D,z(Op,/Oy) 4) it will b 1 ; 1 in further iderati — 
#4 it will be neglected in further considerations. 
of M, = D,(dy,/dy) + D,,(O¥- 2») repeal. - Ses BY faa " —— 
Substituting the expressions (7) into Eqs. (5), the 
he My, = Hy |[(Op,/dx) + (O¥2/2Y)] eee ny ee ee 
" following three equations for the components of defor- 
3 eae mation are obtained: 
~ Oy —< Wy = Ft = Ot - ow oO” 
D, — + D,, A, a ~ alee A) oe 
ox Oxdy Oxoy Oy? Ox ot? 
ic, = Oy ~ Ov 5 OV: 7, OW 5 Ow o”y 
D, — + D.,, —— + An + Ay =, — Ry (w+ = [— (9) 
n Oy" Oxdy Oxdy Ox? Ov Of? | 
_ Oy _ OY, _ Ow » ow "Ww 
ea t+ Kh, 4+ Kh, — +8, —4+¢=-"Uu— 
Ox Ov Ox" Oy? or” 
733 
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The coefficients D,, .. 
(9) determine a unique solution if on the boundary one quantity of each of the pairs 


(Mi, i), (Mi, ¥), (Qh w); 1=%,Y J=%*%,Y 
are given. For the simply supported plate, the following boundary conditions will be assumed: 
or x = 0, M,=0, %|¥=0, w=O0) 
for y = 0,8, M,=0, ¥,=0, w= Of 


Determination of Dynamic Response 


., A, are now complex, as they contain the complex moduli of elasticity in the form (4). 


Eqs. 


In order to determine the response of the plate to a loading random in time and with known space distribution 


g(x, y) it is necessary to determine first the response of the plate to the loading 
° = twl 
q(x, ¥, t) = g(x, ye 


Employing Fourier-series expansion, the loading g(x, y, ¢) may be represented in the form 


n= x u=« 
: : wl 
q(x, ¥, t) = = +S Ona SI Ap, % 9119, e 
m=1 n=1 
where Qn = Mr/a : B, = nx/b 
oe oa: e= 1,2 S,... 


with the coefficients Q,,,, properly determined for each particular type of loading g(x, y). 


of the unit pressure uniformly distributed over the plate, i.e., g(x, y) = const. = 1, then 
Qnn = 16/mnx?; C= 1 SB ae a ao 
The solutions for y,, y,, w are assumed in the forms 
iy = = > u Vimn COS Am Xx sin B,y e” = yp, e™ 
m 
% =D > Qumn SiN am x cos B,y e™” = zy, e™” 


m n 


w = d dX Wan sin an x sin Bye = @ ec 
n 


m 


which satisfy the boundary conditions (10). Substituting (14) into Eqs. 


equations is obtained for Vimar, Qynny Winn: 


(Dan? + FB? -+- K, 
(Dy mB n =) © TT, 0mBn)V 
K ,0nVrmn + RBiVymn = | K itn i KB. 


ined Tw) Vimn + (DryamBn + H, yOmBn)V ymn + K, am mn ~~ a, 
xmn ig (DyBn* - TT ,y0¢m? + K yo Tw? )Vymn + K ‘SV mn = 
— Mw?)Wmn = QOnn 


(11) 


(12) 


If the loading is in the form 


(14) 


(9) the following system of linear algebraic 


(15) 


The coefficients of these equations are complex, and the solutions are, in general, also complex functions of w or jw: 


WV emn(1w), Vymn(tw), Wnn(iw). It is easy to check that the solution for VW,,,,, is 
Elw? + F 


1 
Wan = —_ g — M ” 
. FE + @/w? + ® TS " | 


The functions V,,,, and Vy», can be expressed in terms of W,,,, in the following way: 
K,a,Jw? + @ 
I?wt + Qlw? + B 


KB, 1w? D 
hig Ri ee 
I?wi + Qlw? + B 


V, mn = ] J a n 


The notations in the expressions (16), (17), and (18) are: 


Q@= —( Dein’ + DB? + Tryon? + H,,B,? + K, + K,) | 
@® = (D hm? + A,B x? + kK )(D, Wn? + TT py 0m? fe R) _ (D,., + Hi. Nin bet 
e = —R.a(D,8,! + H,0,? + R,) + R,Ba(DyycimBn + HeyOmBn) | 
Dp = — KB, (Dam? of HTB? + R.) + RK 0m(DyamB T IT, mB) 
& = &K.*a,° + K, “2, 
§ = —H,,Kfan' — H,,K,B,' + (—D,K2 + 2D,,K,K, + 2f,,K,K, + DR?) an?B,? — | 
K2R yan - K,K,8,. 
GS = Kap? P+ K,B,? 


(16) 


(17) 


(18) 


(19) 


(20) 
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ORTHOTROPIC SANDWICH PLATES 7 
If the frequency-response functions JY, = ¥,(iw), ¥, = ¥,(iw), ® = M(iw) are known, the frequency-response 
function of any component of deformation or stress can be easily calculated. 
Considering the stress o,, in the facing, 
o,, = E’ h Ov, ze h Wy 9 
ies "2 Ox "2 dy aaa 
Consequently, the frequency-response function for ¢,, is 
a & ' ; a @ ; 
o,-(to) = — E’., m= > >. Veuntn Sin a, x sin By — E's, e YDS VymrBn SiN Am X SiN BrY = 
ad m w = m n 
lies - ane R: re 
— Dd | 2’ (1 + tn’rz) 5 Vemn@m + FE’ + tn'sy) 5 VumnBn | Si a X sin Byy (22) 
m nu = _ 
and 
a oat ss = ~ os 5 
max Grr(iw) = — DY D>) | E’re(1 + t's) 5 Vamn@m + Ey + t9'ry) 5 VumnBn (23) 
m nu - —_— 
To determine the frequency-response function of the shear stress in the core 6,,(/w), (or: = G,-(1w je"), it is necessary 
to start from the relation’: 
or, = G,[vz + (Ow/Ox) | (24) 
which implies 
&,(iw) = G,(1 + itr) >, Yo (Vann + WirnQ@m) COS Am xX sin By (25) 
m ” 
and 
max 6,,(iw) = GCL + iy) > (Vann + Winnaem) (26) 
m n 
Similarly 
: a EF! a h y kk’ in! h W i - SI ! 27) 
6,(w) = — = > Ly oyll + 2 yy) . ViymnBn + E'r(1 + 19’ xy) = VemnQ@m | SIN Am X Sin By (27 
m n aad —_ 
and Tys(tw) = G,(1 + ity) >, 2 (Vymn + WinnBn) SiN Om X COS Bry (28) 


m n 


Example 
For the numerical evaluation of some of the above 
equations, a plate with the following dimensions has 


been assumed: 
a = 20h, 6 = 16h, f = 0.05h 


The elastic properties of the materials of the core and of 
the facings are described by the set of elastic constants 


hn = f= 

Ey, = Ey = 0.3E 
CoS fo Se = ee 
I 


ie = Pe = 0.3E’ 
G, = G, = 0.3846E 


The shear damping of only the core has been considered, 
with the damping coefficients 


Me = My = w = OI 


Fig. 3 shows the response of this plate to a loading 
uniformly distributed on the surface and harmonically 
varying in time with the angular frequency w. The 
absolute value of the dimensionless deflection ®@D,/gh* 
is given as a function of the dimensionless frequency 
w/wo, Where aw? = D,/h*J. Some insight into the dy- 
namic behavior of this plate may be obtained from Fig. 





!, which represents the deflection @ caused by the 
loading distributed sinusoidally over the surface and 
harmonic in time (in plotting the diagram, the effect 
of damping was neglected). As should be expected, 
three “‘resonant”’ frequencies result from the theory em- 
ployed for this calculation. The classical plate theory 
would give only one of these frequencies. 


Concluding Remarks 


The expressions derived for the frequency-response 
functions of deformations and of stresses may be used 
to determine the behavior of a plate under normal 
pressure which is uniformly distributed over the surface 
of the plate and random in time. In particular, it is 
relatively easy to determine the power spectra of the 
components of deformation and of stresses for an arbi- 
trary power spectrum of loading, even if this is not 
given in the form of a mathematical expression (it may 
be given, for instance, in the form of a diagram). 

The material damping is considered in a manner 
which seems to be as exact as it can be without the use 


of nonlinear equations. The analysis neglects the inter- 


(Continued on page 752) 
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Magnetohydrodynamic Free-Convection 
Pipe Flow 


Kenneth R. Cramer 
Thermomechanics Laboratory, Aeronautical Research Laboratories, 
Air Force Research Division, Wright-Patterson Air Force Base, Ohio 


December 2, 1960 


Tl HAS BEEN SHOWN! that transverse magnetic fields of practica 
strengths exert considerable influence on liquid-metal, free- 
convection, vertical, flat-plate and parallel-plate flow fields 
The extent of influence was determined by the magnitude of a 
nondimensional parameter A which is the ratio of the Hartmann 
number to the fourth root of the Grashof number, 
measure of the relative influence of the magnetic 
In this note the steady, fully developed, 
a fluid of electrical conductivity 


and is a 
and buoyant 
forces laminar, free- 
convection flow of oy through 
a fully 
pipe located in a 
analyzed in terms of the same parameter. 
its influence on the velocity and temperature profiles, the surface 
and the volumetric flow rate is deter- 


constant-temperature, vertical 
field of strength Bo is 


The magnitude of 


submerged, open-ended, 


transverse magnetic 


shear and heat transfer, 
mined. 


ANALYSIS 
The magnetohydrodynamic, free-convection, — pipe-flow, 


boundary-layer momentum and equations and_ their 
boundary conditions at the inner pipe radius a |Eqs. (1) to (3)]} 
are transformed to the nondimensional equations (5) through 
(7) by the introduction of the new variables, Eq. (4), where NV 


is the pipe temperature per foot and Z is the pipe length. 


energy 





dy L du a By? rT ~ T 
dr? r dr pv line ae (] 
(d?T/dr?) + (1/r)(dT/dr) = O (2) 
r=a:u=O0O8&T = T°F (3) 

u = (gBN)'?Lf(n), n = (gBN)'4yp-V2y ) 
a = pl2(gBN)- "4, e = (7 —TINT. =Ts (4 

Tw — To = NL, A? = a By?/p(gBN)'? \ 
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NONDIMENSIONAL PIPE RADIUS |7| 


Fic. 1 Velocity profiles 
id + (] nf’ — Atf = —@ (5 
0” + (1/n)0’ = 0 (6) 
n = 1:f(1 0 & a1) = 1 (7 


The solutions of Eqs. (5) and (6), where the velocity is in terms 
of modified Bessel functions of the first kind, are then 


¢6= 1 (8 
f 1/A2}1 — [Ip(An)]/[o(A)] 5 (9) 


A separate solution, Eq. (10), is required for vanishing A since 


it appears in the denominator of Eq. (9 
fio = ~(1 — 2? (10 
A= } 

The relations [Eq. (12)] for the volumetric flow rate Q are then 

obtained from Eq. (11) after substitution of the appropriate 


velocity functions: 


(12 


RESULTS 

The solution of the energy equation is independent of the 
influence parameter; thus the temperature field is not altered 
by the magnetic field. In addition, no heat is transferred be- 
tween the wall and fluid in the fully developed, steady-fiow, 
constant-wall-temperature case. 

Figs. (1) and (2) demonstrate the influence of the magnetic 
field on the velocity profile, surface shear, and flow rate. In- 
creasing the influence parameter produces the expected flattening 
of the velocity profile and corresponding reduction of the surface 
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Fic. 2. Surface shear and flow rate 


Their asymptotic trend for large values 
separation and flow stoppage 


shear and flow rate 
of A demonstrates that flow 
cannot be achieved, since the magnetic field exerts only a re 
tarding force on the flow field 

Magnetic fields of practical strengths can exert considerable 
influence on free-convection pipe flows. For example, from the 
definition of A, a magnetic field of 10,000 gauss will produce an 
influence of A equal to ten on the flow of mercury through a 1-ft 
pipe with a wall temperature of 100°. Less than 1,000 gauss 
would be required to produce the same influence on liquid 


sodium and potassium flows 


REFERENCE 
Several Magnetohydrodynamic Free-Convection Solutions, 


1 Cramer, K. R., 
to be published 


A Mass-Transfer Finite-Difference 
Formulation Employing Crocco Variables* 


J. P. Moran and P. B. Scott 


Naval Supersonic Laboratory, Massachusetts Institute of Technology, 
Cambridge, Mass. 
December 8, 1960 


fp INJECTION of a foreign gas through a porous wall into the 
boundary layer is being studied as a possible means of 
cooling a high-speed vehicle. The complexity and nonlinearity 
of the describing equations makes it impossible to solve them 
exactly except for cases when certain restrictive boundary condi 
tions apply. Here the laminar, binary, boundary-layer equa 
tions for the case of helium injection as derived by Baron! have 
been solved by a finite-difference scheme allowing for more 
general boundary conditions. This formulation employs the 
Crocco form of the equations using the same finite-difference 
approximations as Fliigge-Lotz and Baxter.2 The descriptive 
equations were programed for solution on an IBM-704 computer 
permitting wall temperature, injection rate, and external-flow 
properties to be varied with distance along the body A finite 
difference formulation of this type is not adequate when velocity 
overshoot exists—i.e., somewhere within the boundary layer 
the velocity exceeds its local external value—-since the dependent 
variables have infinite derivatives at the point of maximum 
velocity and are double-valued in the vicinity of overshoot 
For example, helium-injection boundary layers at moderate 
levels of injection are free from velocity overshoot only when 
the external pressure gradient is at most only slightly negative 
This restriction therefore must be placed upon the external- 
velocity distribution in that case 

Solutions have been obtained for the problem of a porous 
plate having a solid leading edge, the downstream porous region 
having a mass-injection distribution varying inversely with the 
square root of distance from the tip. The numerical analysis 
carried out in the detailed report® gives reason to conclude that 
such solutions are a quite accurate representation of the exact 
solutions but for a small region immediately downstream of the 
solid edge. In Fig. 1 the variations of skin friction, heat transfer, 
and surface concentration are presented in normalized form 
against distance from the start of injection. The parameters 


employed are: 


crv Re, StxV Re, and 1, 


in which c; = skin-friction coefficient, Re Reynolds number 
based on stream properties and length from tip, St, = Stanton 
number, and ¢;, = mass concentration of helium at the surface 


The distance is normalized with respect to solid-leading-edge 


* This research was supported by the USAF under Contract AF49(638)- 
245 monitored by the Air Force Office of Scientific Research of the ARDC 
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Fic. 1. Normalized skin friction, heat transfer, and surface 
concentration vs plate distance 
length—i.e., ¢ = x/L. The normalized forms are: 


/ 
(c¢V Re)g=1 — cyV Re 


pr 7 j 
(c;V Re)g=1 — (c7V Re)z20 
(St.V ‘Re)e=1 — StxV Re 
or = ] ] 
(SthkV Re)ga1 — (StaV Re)zug 
Clu 
oe = 
(Ciu R= a 


The asymptotic values for the parameters were obtained from 
the similarity solutions of ref. 1. Boundary conditions for the 


solutions are given in the following table: 


( pv )u Tu F —> 
on / Re : = fw for LO 
Run p'u M T stream (vet. 1) 
1 0.1083 0 0.8 —0.193 
Z 0.1597 0 0.8 —0.300 
3 0.1986 0 0.8 —0.373 


fw is the injection parameter! and is related closely to 
[( pv) w]/( pu") V Re 


A solution for the boundary conditions of Run 3, but with the 
inclusion of an external Mach number of 3.25, was also obtained. 
The resulting wall concentration values were very close to those 
of Run 3. Shear and heat-transfer results were significantly 
different from Run 3 but are not presented here since exact 
asymptotic values are not known. 

The curves are applicable to a solid-tipped cone after con- 
tracting the dimension ¥ by the 1/3 power in accordance with 
Mangler’s transformation. A cone of this type is useful in 
obtaining an experimental check for the similarity solutions of 
ref. 1 and others provided the solid tip is sufficiently small so 
that similarity is attained somewhere along the cone. The 
cone tip ideally should be porous, but structural limitations 
require that it be solid. The solution presented in Fig. 1 shows 
the rate of decay of the solid-leading-edge influence over a portion 
of the region of influence. The solutions were not carried to their 
asymptotes due to the prohibitive computer time required. 
However a combined application of the finite-difference program 
and Covert’s approximations‘ as performed in the full report 
indicated that similarity profiles are approached along a cone 
within three solid-tip lengths downstream from the start of 
injection, 
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Heat Transfer in the Incompressible Boundary 
Layer on a Flat Plate With Arbitrary Heat Flux; 


A. G. Smith* and V. L. Shah** 

Department of Aircraft Propulsion, The College of Aeronautics, 
Cranfield, England 

December 1960 


SYMBOLS 


distance from the leading edge 
= distance from the leading edge where heating starts 


Ml 


a 

é = variable of integration 

p = density of the fluid 

Cp = specific heat 

K = thermal conductivity 

v = kinematic viscosity 

h = coefficient of heat transfer 

U = velocity in the boundary layer 

U: = mainstream velocity 

t = temperature in the boundary layer 

h = mainstream temperature 

tu = wall temperature 

6 = hydrodynamic boundary-layer thickness 
A = thermal boundary-layer thickness 

dw” = heat flux at the wall (quantity of heat per unit area and time) 
m = slope of ramp heat input 

Cy = friction factor 

Pr = Prandtl number 

Re = Reynolds number, Uix/v 

Vu = Nusselt number, hx/K 


ECENTLY Rubesin,* Eckert,? and Reynolds ef al.' have derived 
analytical equations to find the heat-transfer coefficient on 
a flat plate for an arbitrary varying wall temperature. Following 
Klein and Tribus,‘ they have derived equations to find wall tem- 
perature for an arbitrary heat flux, from the heat-transfer equa- 
tions for an arbitrary wall temperature. 
Equations derived below directly from an integral energy equa- 
tion are simple and still give results within 1.5 percent of those of 
Reynolds et a/.! and Rubesin.* 


Laminar Boundary Layer 

For a step heat flux at x = a, such that g.” = 0 for x < a and 
Gw” = b for x >a, the total energy at any section x(x > @) is given 
by 


A U 
Gw"(x — a) = pU; Cy(tw — t)6 f, (4 Jaa — 6)d(y/6) (1) 
1 


For assumed cubic velocity and temperature profiles, 
(U/U1) = 3/2 y/6) — 1/2(y/6)8 (2) 
0 = (t — tw)/(t — tw) = 3/ey/A) — 1/e(y/A)3 (3) 


Eq. (1) becomes 


- i 3 /s\? 3 (a\ 
Gw'(x — a) = pU,Cj(tw — t1)6 20 \3 — 380 \3 (4) 


We are treating the case where the thermal boundary layer is 
thinner than the velocity boundary layer—that is, (4/6) < 1. 
The second term in the right-hand expression of Eq. (4) is then 
small compared with the first, and can be neglected. The heat 


transfer at the wall is given by 


inn” = Wty — t «(2 °e=e—= w 
Gop = Uly — == = (2 
” dy),-9 2 A , 


T This is a part of the research supported by U.S. Air Force under Contract 
No. AF 61(052)—267. 

* Professor and Head of the Mechanical Engineering Department, Uni- 
versity of Nottingham, Nottingham. 

** Research Assistant 
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Fic. 1. Step and delayed-ramp heat input (laminar boundary 
layer) 


For a laminar boundary layer, the velocity boundary-layer thick- 
ness is given by 
5 = 4.64x/V/ Re (6) 
Thus for a step heat flux, from Eqs. (4), (5), and (6 
2.395b( Re)?! Pr)?/3 | 
(ty -— hh) = rr il — (a/x)j 
plc p 
and Nuz/V Re = 0.418(Pr)/3[1 — (a/x)| 7 (8) 
Here 6 is the magnitude of the step heat flux. When compared to 
the similar equations of Eckert? and Rubesin,* it can be seen that 
Eq. (7) is relatively simple, for it does not contain the incomplete 
beta function. For any arbitrary heat flux, the wall temperature 
and the heat-transfer coefficient can be calculated by superim- 
posing small steps. Thus for any arbitrary heat flux 


2.395( Re )!!2( Pr)2/3 x dqu "(£) 
/ ¢ 1 wae rt 
( —-i)= (1 — &/x)!/8 dé (9) 
0 ¢ 


f 
pUiC, dé 
and 
Nuz 0.418( Pr)*/3g,."(x 
= = (10) 
x 
V Re f (1 — &/x)!!3[dqw"(&)/dé]dé 
0 


It should be noted that the integral of Eqs. (9) and (10) must 
be taken in the Stieltjes sense. For a delayed-ramp heat input, 
when gw” = 0 for x < a, and gw” = m(x — a) for x >a, Eqs 


(9) and (10) give 


1.797( Re)!/%( Pr)?!3mx a 
(te — hh) = - {1 — a/x]*'3 (11) 
pl i¢ 
and 
Nuz/~/ Re = 0.557(Pr)!/3[1 — a/x] 71/3 (12) 


Turbulent Boundary Layer 
For a step heat flux in the turbulent boundary layer, we as 


sume velocity and temperature profiles 


(U/U,) = (y/6)!" (13) 
and 
0 = (t — tw)/(tr — te) = (y/A)!? (14) 
Then Eg. (1) becomes 
Gw"(x — a) = pUC,(tw — t)6[(7/72)( 4/6)8/7] (15) 


The heat-transfer momentum analogy derived by Reynolds et al.! 


for step heat flux is 


Ste = Gu" /[(tw — tr)pUsCp] = (Cy/2) 4/8) 17 (16) 
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where the Prandtl-number dependence factor is applied after 


ward. 
For a turbulent boundary layer 


C; /2 = 0.0296( Re) ~°-2 (17) 
6/x = 0.37(Re)~°%? (18 


Thus for a step heat flux, from Eqs. (15) to (18), 


32.95b 


(te —~—h)= —— (Re)®-% Pr)®-*{1 — a/x]!"9 (19) 
pC, U) 
and 
Nu,/(Re)®’ = 0.03035( Pr)®-*[1 — a/x] 7/9 (20) 


Here 0 is the magnitude of the step heat flux. A similar equation 
of Reynolds et al.' for a step heat flux also contains the incom 
plete beta function. For any arbitrary heat flux, the wall tem 


perature and the heat-transfer coefficient are 


32.95/( Re)®-*( Pr)®-4 x dgjw"(E) 
(tw — 4) = ane (i = €/a;""" dé (21) 
pCpU, 0 dt 


and 


Nu, 0.03035( Pr )°-®Gw" (x) 


(Re)®-8 x dgjiw"(&) 
(1 — &/x)*'® dé 
0 dé 


In Eqs. (21) and (22) the integral must be taken in the Stieltjes 
sense. 

Results of the above equations are compared with those of 
Reynolds et al.,! Eckert,? and Rubesin,’ for both turbulent and 
laminar boundacy layers, and are plotted in Figs. | and 2. For the 
laminar case, the results of the above equations agree very well 
with those of Rubesin,* but not with those of Eckert.2 The dis 
agreement with Eckert is merely a numerical factor. Eckert’s 
(tw — t,) values are 8 percent lower than those of Eq. (7). For the 
turbulent case, the results agree within 1.5 percent of those of 
Reynolds et al.! 
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Beam Column of Nonuniform Sections 
by Matrix Methods 


H. Saunders* 

Missile and Space Vehicle Dept., General Electric Co., Philadel- 
phia, Pa. 

December 10, 1960 


SUMMARY 
Aun analytical solution of a beam column of nonuniform cross 
Most of the 
available analyses consider a varying or uniform section but very 


section was conducted and adapted to matrix form. 


few take into account abrupt changes. The matrix procedure as 


applied herein is called the ‘‘transfer-matrices’’ approach. ! 


NOMENCLATURE 
(EI); = bending rigidity at station 7 
Vi = deflection at station i 
P = axial force 
Mi, = moment at station 7 
0; = shear at station 7 
x = distance along beam 
Leas = distance along beam between stations i andi + 1 
FY = Laplace transform symbol 
F; = lumped static load at station 
ANALYSIS 


A“ ELEMENT of a beam column with moments, shears, and 
axial forces acting upon it is illustrated in Fig. 1. The 
beam is divided into finite discrete stations with the static loads 
lumped at the stations. The differential equation for a typical 
ith section is 


d?y 
(ED in — = —Piy — vi) + Mi + Qi(x) (1) 
dx? 
> 
Simplifying, dividing by (£/)j+; and letting E = &+,%, one 
(EI) ;+1 
finds 
COs & piliet — &;,; sin @i 11; 
sin @; 10; 41 
cos @isalins 
@ix 
Aisa = | 1 = cos @ighsas sin @ailiss 
@i12(El) 541 ain(El); 41 
@igiliga — Sin @aylicy 1 — cos @yilias 


&i41? @.17(El) 541 
The following relations can be written at station 7 + 1 between 
(¢ + 1)" and (i + 1)4 


yitr® = yor” 


/ > / 
yaty® = yey! 


and similarly for stations i + 2,7 + 3...n. 
in terms of station 1, 


Q Ay Lie 113 0 Q 
M\ _ |Aa Ax fo; 0 M 
y'|  ~—s | As Ax A33 0 y’ 
y \n Ay 1 43 Ay 1 y | 


For the various necessaty conditions of the beam column: 


(a) Fixed ends: yn, = yw = yn’ = m1’ = O 
Determine the critical load from A3;Ay. — AyAz = 0. 
(b) Pin-ended: y, = M, = y, = M, = 0 

The critical load is AnAqg — Ay»Agy = 0. 


(c) Cantilever: yn = yn’ = Qi = M, = 0 
Obtain the critical load from A;3; = 0. 
(d) Clamped at station n and pinned at station 1: 
1 = M, = yn = yn’ = O 
* Engineer. Also, Lecturer, Dept. of Mechanical Engineering, Drexel 
Institute of Technology, Philadelphia 
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d*y ie ee M; Qx 

itty = B49; t (2 
nwo ee ee 


Applying operational methods by taking the Laplace transform? 
of Eq. (2), then, one obtains 


in P i aAi+ i 
(s? + ai+12)y(s) = sy(O) + yO) + + 7 


s*E] iH 


The initial condition are at x = 0, thus y(0) and y’(O) are actually 
y, and y;’, respectively. 


Dividing by (s? + @;+,?) and taking the inverse, 


Mi Sin &i4ix Mi(1 — cos &+1x 
y(x) = yi + ae 
Qi+1 @j+17(E1) i+) 
QOi(@itix — sin &+1x 
= (4 
@i+1° EI it+1 
Differentiating Eq. (4), 
; . M; sin @+ix | Qi(1 COS &i+14 : 
y' (x) = yi COS Bix + i t rye (5 
@i4+1(E/) +1 &+,?(E1)i+1 
Differentiate Eq. (5) and divide by E/ for the moment. Simil- 


arly, differentiate the moment equation and obtain the shear. 

Now consider that x = /;+; and thus y, y’, J and Q immediately 
to the right of station 7 + 1 can be expressed in terms of the cor- 
responding quantities at station 7. The resulting deflections, 
slopes, moment and shear are expressed in matrix form: 


y| atiy® = | Asi] ly i (6) 
where 
yl ei® = Q and iy ; ¢) 
\ M 
y’ y’ 
y oh y 


and A;4; is the transfer matrix.! 


+1 —@: 1 2(El);4; cos &yiliss 0 
— &;,;(El);.; sin @iailias 0 

(7 
cos @issligs 0 


— sin @ailias 


Ai+l 





es A ae 
(e) Pinned at station n and clamped at station 1: 


The critical load is obtained from A 3)A. 43 


yn = M=n=n =0 
AwyAy = 0. 


The simplest approach to the matrix solution is to assume 


Determine the critical load from 4,Ay. — 








Y | 
re Yiti Yi 





x | 
Fic. 1. Element of beam column. 
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> = 


ai+i1; then ai+e can be selected from @+;2 = &+2?{(E/)i+2/- 


El i+]. 
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Further Note on Airfoil Theory in 
Hydromagnetics* 


G. S. S. Ludford 

Div. of Applied Mathematics, 
Brown University, Providence, R.! 
December 19, 1960 


HE PREVIOUS NOTE! was concerned with the lift loading on a 
pee airfoil moving at right angles to an applied magnetic 
field in a fluid of infinite electrical conductivity. The object of 
this sequel is to point out that the three basic integrals (J), Jz, Js 
below), which were there expressed in terms of hypergeometric 
functions, are in fact simple algebraic functions; and to show 
that this leads to results comparable in simplicity to those in 


ordinary airfoil theory 


The airfoil has camber line y = V(x), —b < x < 6, and moves 
in the negative x-direction with velocity U, the magnetic Mach 
number being m = tan Bz, 0 < 8B < 1/2. The corresponding 
results in ordinary airfoil theory are obtained by setting m = —, 


B 1/2 in the following 
First note that 


° ( +s\8 ds 
~» Vn t-s 
T b+ x) 8 
. 1 — cos Br ) 
sin Br b-—x 
Thus the lift loading! 


271+ m2? sin Br [b — x\8 
h(x) = > pU? ( I(x) + cos Br 
m 7 b+x 


on a flat plate at unit incidence [¥’(x) = —1] becomes 


L(x) = 


I(x) = 2 cosecBr pU2[(b — x)/(b + x)]F 


The total lift on the plate is 
b 
Io = l(x)dx = 48x cosec? Br-pU*b 
—b 


while the pitching moment about the center is 
* This research has been supported by the Office of Naval Research, under 


Contract Nonr 562(07) with Brown University 
T This can be seen by transforming the argument of the hypergeometric 


functions to its reciprocal As a by-product we obtain the complex velocity! 
0¢/0x — 10g/dy = il’ [(s — b/s + b)3 — 1] for a flat plate at unit incidence. 
r y 
P 
J 
/ 
J 
i 
4 
y, ‘ 
9 F Fic. 1. 
. 5 
8 J 
Ys 
a 
ri 
27 / 
6 | 
352} 
429 3 
- 
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xl (x )dx 


=" 


43? x cosec? Or: pl 242 


The lift force acts at x = —Sb—i.e., at the 8/2-chord point for- 
ward of the center; it moves toward the center as m decreases 
i.e., as the magnetic field becomes stronge1 

The figure shows a graph of the lift coefficient C,, = Lo/pl%b 
482 cosec? Bx versus 1/m. As m decreases from its value © in 
the nonmagnetic case, C,, first decreases from 27 6.28 to a 
2.33 (1/m = 0.429 


to 27 again for m = 1 and indefinitely thereafter as m tends to 


minimum of 5.52 for m It then incieases, 


zero. 
Secondly, note that 


b—x\? dx 
To ‘y= - 
=" e+? Ss 5-23 


x h «\8 
—_- 28b — § + COS Br ( § 
sin Br O+S 


Thus the total lift and the pitching moment (about the center), 


"gee , b 
1 + m? : sin Br b+s\? : 
L= _»* put (| ) In(s) ¥'(s)ds 4 
h 


m1 Oo-S 


1 + m? sin Br {°° b + 5\8 
wet pU? ( In(s)¥"(s)ds 4 
ml T , b-—s 
eh 
cos Br | v V'(x ax| 
r= 


experienced by a general airfoil, become 


b b+s\8 
L = —2 cosec8z pl? ( ) Y’(s)ds ML 
_p \O-s 
and 
b b + s5\8 : 
M = —2 cosec8x pU? (28b — s)¥'(s)ds = \2Mo 
-. ew se 


When the same airfoil is set at an angle of incidence a these 
values are augmented by aly and aM respectively. The 
moment about the 8/2-chord point is then independent of a and 
this point is the focus? of the airfoil. For a = —d, there is no 
lift, while for a = —d, there is no moment about the center 
Thus \, and 22 define the first and second axes.?, Expressions for 
them directly in terms of the shape of the camber line are ob- 
tained on integration by parts 


sin Br b V(sjds 


A = : 
: (b — s)!t8(p “—P 


. —b  * 


6 (482 — 1)b2 —2 Bbs + s? 


ie 7 287b? 


We have used the fact that V(x) vanishes like ) — « at the trailing 
edge x = Bb. 
In terms of the modified Glauert series,' given by 


— U(tan 6/2)! ~ 28 ¥’(b cos 6) = Bo + } By, cos né, 
n l 
x = beos #0 LS A< 7) 


they become 
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4, = (sin Br/4BU)(2By + B,) 
A» = (sin Br/8B2U)(2(48B — 1)Bo + 2(28 — 1)B, — By] 


These expressions contain only the first three Fourier coef- 
ficients, and the remaining coefficients therefore correspond to 
changes in the shape of the camber line which have no effect on 
the lift or moment. Of course the results given above may also 
be derived by means of the modified series, in the same way that 


Glauert used his original series.* 
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Determination of Real-Gas Stagnation 
Temperature Based on Mass-Flow 
Consideration 


Frank L. Clark and Charles B. Johnson* 
Langley Research Center, NASA, Langley Field, Va. 
December 21, 1969 


O~ OF THE MAJOR PROBLEMS encountered in the operation of 
high-enthalpy test facilities is the determination of high 
stagnation temperatures (J > 4,000°R). Since these tempera- 
tures normally exceed the capabilities of available thermocouples, 
it is often necessary to resort to one of the various optical meth- 
ods—such as pyrometry or the sodium-line reversal technique 
for the determination of these temperatures. A method has been 
presented! whereby real-gas stagnation properties in a_high- 
enthalpy system can be determined if the mass-flow rate through 
the system and the stagnation pressure are known. In the present 
analysis, this method was applied to determine stagnation tem- 
perature for stagnation pressures ranging from 100 to 1,000 atm. 
Nitrogen is considered to be the flow medium and the expansion 
of the gas from the stagnation chamber is assumed to be an 
isentropic and steady-state process. The maximum tempera- 
ture dealt with (5,000°R) was low enough that the gas remained 
in an undissociated state; however, the method may be extended 
to other gases which may or may not be dissociated, provided the 
thermodynamic properties of the gas in the desired pressure and 
temperature range are known. 

A one-dimensional analysis of the flow is used to evaluate the 
mass flow per unit area at the throat of a nozzle for various stagna- 
tion pressures and temperatures. The evaluation of mass flow 
per unit area at the throat of a nozzle is a trial-and-error pro- 
cedure and was performed in the following manner. Starting 
with an assumed stagnation pressure and temperature, an initial 
value of stagnation entropy was obtained from published data.?: 4 
The flow in the system was expanded by reducing the stagnation 
temperature at a constant value of stagnation entropy. At each 
assumed temperature the values of pressure, density, and en- 
thalpy were obtained from Ref. 2 or 3; the velocity at each point 
was calculated from the energy equation. 

The result of the analysis is presented in Fig. 1, wherein mass 
flow per unit throat area is plotted against stagnation tempera- 
ture. The dashed lines represent ideal-gas calculations. The 
spread between the ideal-gas and real-gas curves becomes larger 
as the stagnation pressure is increased. It may be seen from 
Fig. 1 that a significant error in stagnation temperature could re- 
sult by using the ideal-gas relationship. For example, at a mass 
flow per unit throat area of 16,000 Ib/ft?-sec and a stagnation 


* Aeronautical Research Engineers 
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Fic. 1. Variation of mass flow per unit throat with stag 
nation temperature for various stagnation pressures 


pressure of 1,000 atm., the values of stagnation temperature as 
given by the two curves differ by approximately 600°. 
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Drag Reduction of an Airfoil by Injection of 
Sound Energy 


Paul K. Chang 
Associate Professor, Dept. of Mechanical and Aeronautica 
Engineering, The Catholic University of America, Washington 17, 


December 19, 1960 


HIS PAPER presents preliminary results on the reduction of 
total drag of an airfoil by injecting sound energy into the free 
airstream. It was found that the maximum effect was assured 
when the sound energy was directed opposite to the free stream 
and it is shown that the drag can be reduced by more than 2() 
percent at zero angle of attack. At an angle of attack of 15 
the power saved by the drag reduction is 19 times as great as the 
power input required to generate the sound energy. 
Measurements of total drag with and without sound energy for 
a wooden wing of rectangular planform, aspect ratio of 6, chord 
length of 3 in., and maximum thickness of 15 percent at 0.6 chord 
length were made at the Catholic University of America.' The 
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TABLE 1 
Ratio of Power Reduction to Power Input 
; Velocity of airstream, Frequency of sound, Strouhal Power 
\ngle of attack Total Cp reduction ft/sec cps number reduction/power input 
0° 0.007 35.7 1,500 17.52 6.7 
15° 0.022 34.5 2,000 24.20 19.0 
range of Reynolds numbers based on chord length and free-stream 2 
velocity was 8.34-8.7 X 10‘, the range of frequency was 500- 3 
: e: % ry 
2,000 eps, and angles of attack varied from —15° to +15°. a 
Sound energy from a Hewlett-Packard audio-oscillator was in- 3 
iected into the air stream in the range of Strouhal numbers from 3 
5.82 to 24.2. The location of sound source with respect to air s ae 6° 
. . . . . e 
stream is important, and in order to obtain the maximum effect fo) nome a 
‘ ; , f ete +15° 
of sound, the audio-oscillator was placed behind the airfoil, so z ° =e 
that the sound was propagated against stream direction 3 
Some of the test results are shown in Figs. 1 and 2, the data od 
being average values of several test runs. The curves show that 5 "4 
the drag at various angles of attack can be reduced over a range o o 
of Strouhal numbers. Strouhal number is defined as S = nc/V, " 
where 7 is the frequency of sound (cps), ¢ is the chord length of an ° ; 
airfoil (ft), and V is the free-stream velocity (ft/sec). It should =z o a a 
; a d xz wf asese 
be noted that in these preliminary results the Strouhal number o 
was varied by changing only the frequencies from 500 to 1,500 eps, g aeo 
but the flow velocity was 35.7 ft/sec for all frequencies of air flow 3 
tested. The minimum total drag was measured in the range of “4 4 
S = 11-18 at various angles of attack, and the maximum reduc- =z oOo 
tion, amounting to 21 percent of total drag, was attained at S = re) 
a f : - 
17.5 and at zero angle of attack. 
.@) a 8 12 16 20 
ag 
g STROUHAL NUMBER =S=-V& 
@ 
“ ons OF ae Fic. 2. Measured total drag ratio of airfoil vs. Strouhal number 
° S = STROUHAL NUMBER 
Se &@s 
TL = FREQUENCY OF SOUND, 6YOCLES/ SEC. The ratios of power reduction to power input to the audio- 
) A C * CHORD LENGTH, FT. oscillator are shown in Table 1. 
of V © FREE AIR STREAM VELOCITY, FT./SEC. eo 
igh REFERENCE 
Sei 1 Betz, A. J., Jr., Burke, P. J., Hendricks, C. M., Lawrence, J c?., and 
Rogers, L. W., An Experimental Investigation of the Effect of Sound on the 
oO ee Ae-16° Drag Characteristics of a Finite Wing, B.S. Thesis, Catholic Uiv. of America, 
own ——__ 1960. 
ge o oO RS eet 
ech 
5 scence 
Ww © 
o of 
i Comment on ‘‘Shear-Lag Solutions for 
= Sheet-Stringer Panels by Means of a 
° aj=#+ise H d ic A I 59% 
of iar ee ydrodynamic Analogy 
o Oo 
= ‘a a s-10° Martin Goland 
a ee — President, Southwest Research Institute, San Antonio, Texas 
3 4 @o December 24, 1960 
7 < OFT 
; lool ro) 
oO 7 . 
= oe M‘ M. Fine of Dowty Rotol Ltd. has kindly called my atten- 
s= ° ° . . . 
a « - tion toa paper which he coauthored, entitled A Method of Es- 
es 410° : . . + . 
z oo oor timating the Direct Stress Concentration Round Holes in Reinforced 
of o eee ” Aus Be Sheet, by M. Fine and A. Pellew, A.R.C.R. and M. No. 2604 
aed —— aoe (Gt. Britain), May 1942. This earlier paper covers essentially 
" the same material as my own treatment, and my apologies are 
ni due these authors and the readers of the JOURNAL for this dupli- 
»() 4 - : : . abe : : 
; ° cation. My only excuse is the growing difficulty of conducting 
0 4 : 2 16 20 comprehensive searches of all literature sources 
1€ ‘ : 
ne My own version of the problem does include a more careful 
STROUHAL NUMBER S*——<— ewanere se aly: 
Vv treatment of the boundary conditions, permitting a discussion of 
Tr 
d Fic. 1. Measured total drag coefficient of airfoil vs. Strouhal * Journal of the Aero/Space Sciences, Vol. 27, No. 4, pp. 291-295, 303 
d number. April 1960 
e 
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cases when the cutout boundaries are rigidly fixed, as well as 
certain types of cutout reinforcements. In all other respects, 
the two papers show a remarkable similarity of approach. 
EDITOR’S NOTE 

The following communication has been received from Mr. Fine: 
“The authors of R. & M. No. 2604 feel that innocent duplication 
of prior work is unavoidable at present, and are grateful to 
Martin Goland for his acknowledgment of priority. The first 
such use of the Laplace equation, in effect, is possibly by Younger; 
see his Mechanics of Aircraft Structures, 2nd Ed., McGraw-Hill, 


on shear lag.”’ 


Comments on ‘‘On the Rotational Motion 
of a Body Re-Entering the Atmosphere’’ 


Duane F. Dunlap 

Aerodynamicist, Chrysler Corporation, Missile Division, 
Detroit, Mich. 

December 27, 1960 


ETHODS HAVE BEEN DERIVED! to compute the angle-of-attack 
envelope of a body re-entering the atmosphere. Both 
spinning and nonspinning bodies were investigated. Several in- 
consistencies were noted in the derivations and discussion of this 
paper. These will be commented on in the following note. 
In combining the aerodynamic moments with the rotational 
equations of motion, the damping terms Ay and Ay were not 
This results in the im 
The definition should be 


divided by a moment of-inertia term. 
proper definition of the term A y,y. 


Kun = —AmnVeI, (1) 


Without the moment-of-inertia term, the definition of Axsyy 
would be dimensionally incorrect. 

It was also discovered that the constants of integration in- 
volved in the solution of the coupled equations of motion were 


incorrectly derived. The correct constants are as follows: 


3, — BO) ¥ (0) — ay 

a. B,(0) B, = oh )-— a 
I,w,/T: Etts/Ts 

it ay — a,(0) B, = By sad B,(0) 
Tewr : Tw, | P 


A3; = a — afO) — A, B; = Bo — B,(O) — B, 
Ag = B,(O) — Bo — As By = w—- a/(0) — Be 
Another error which seems to compensate for the erroneous 


constants of integration was found in the integration of the damp- 


ing terms associated with the coupled case. The integral 


(1/2)d/dt)V wn? + (Lew, /21.)? dt 
Jo V won? + (Lqnr/21,)? 
when properly carried out leads to the result 

— (1/4) In [w,? + (rer/2I,)] 0 


which must be evaluated for the limits of integration. Since 


w, = Oatt = 0,and w, = w, att = t, weare led to the result 


os f (1/2)(d/dt) Von? + (Lr, /21,)%de, 
0 V oon? + (Leo; /2I2)? 

j (Ie:/2I,)'/? 
In < = Tg 
Veen? + (Lee /212)?|1/4f 
Apparently the proper limits of integration were not applied to 
the coupled case as they were to the planar case. Substitution 
of this term and the A's and B's into the assumed solution will 
yield 
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Aap — (In,\'? 
Aap(O) 7 a) 
e 


[—(m/I,)1zV2e"" (Az + Ap) + (Izws/2I,)!/7}!/4 


—[(A, + Amy)/2k sin et 


This equation is identical to the one published in Ref. 1 except 
for the Ayw term and the negative sign of the exponential. It 
appears, then, that there were compensating errors in the de- 
rivation of the A’s and B’s and in the integration of the damping 
terms. The negative sign which appears in the exponential of 
both the planar and the coupled cases seems to be an oversight, 
otherwise the envelope would diverge. (Both k and sin # are 
negative. ) 

Some comments on parts of the discussion are also in order. 
The statement that A yy is always positive, as well as the impli- 
cation that Cra can be negative, are incorrect. To the best ex- 
perience of aerodynamicists Ayy is often negative and Cra is 
never negative. Negative values of Aywy generally stem from 
pressure effects, rather than purely viscous effects. 

Aside from these comments, the results of the article have beer 
of great utility to the author. 


REFERENCE 
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Author’s Reply 


T. B. Garber 

Electronics Department, The RAND Corporation, 
Santa Monica, Calif. 

January 26, 1961 


N Rer. 1, four specific criticisms are raised in connection with 

the writer’s paper On the Rotational Motion of a Body Re-Enter- 

ing the Atmosphere. These criticisms are discussed below in the 
order in which they appear in Ref. 1. 

1. It is noted that the damping parameters, Ayw and Ky, 
were not divided by the moment of inertia, J,, in Eqs. (46) 
(50).2 
inconsistency. As a consequence the proper definition of A sry, 
Eq. (60) in Ref. 2, is 


As pointed ouvt,! this omission results in a dimensional 


Kyw/I, = —AmnVe* 


Apparently the omission of J, occurred as the result of a notational 
change which was made in writing the paper. The proper value 
of Ayw was used in the preparation of the curves presented in 
Ref, 2. 

2. It is stated! that the constants of integration of Eqs. (74) 
Then it is pointed out that the in 
tegrals of the exponential terms are apparently evaluated im- 
properly? so that a compensating effect is achieved. These 
integrals, as they appear in Eqs. (74) and (75), are indefinite 


and (75)? are incorrect. 


Thus at very high altitudes, 


AMN e* dh 1 / - Izwr\? 
Svat se L-s sin & — 3 In You? + zr; ) 
e e 


»>—(Amve*/k sin &) 


I 


& . 


, Tc. \? | 4 ( 1) 
@n? + 1, 
A Lek “ 1 ° Irwr\? 
_L-s% ® 2 in Vou r oT, ) | 


gn" Axe**/k sin ®) 


: Iewz\?] "4 (2) 
w,? + Py; 


where w; > w, and h = V sin ®. 


S o2 dt 


| 


II? 
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lo illustrate the process involved in determining the constants 
of integration, let aj=;. = ay — as(O) where hy is very large 
e.g., 500,000 ft. All other initial conditions are zero. At very 
large values of h, w, approaches zero, and the exponential terms 
of Eqs. (1) and (2) approach 1, since k is negative. It can be 
shown that, for the assumed conditions, Eq. (74) reduces to the 


following expression 


Tw, \'2 
ag — az(Q) A 7 (3) 


Eq. (3) is consistent with the results of Ref. (2), In Ref. 1, limits 
of integration have been employed in evaluating Eqs. (1) and 


above. Thus 


Ss o; dt Si . . dh 
~~ o ho V sin ® 


é é 
_ AMN (ek! kh \ 
(I,w,/2T,)"%e k sin ® (4) 
(wn? + (rw /2],)?| "4 
: 5 oh (eth — ekh,) 
m dh . ») 1/2 k sin ® - 
Sh V sin @ “" a" (Lrw,s/21,)"* {o) 
[w,? + (eo,/21,)2]" 
The corresponding expression for ag — a,(0) is 
(I tts 27,)} 2 r 

ayo — a;s(O) = 1; A (6) 


(Izw2/2I,)"2 


However, because of the fact that a multiplicative constant is 
introduced when definite integrals are employed, one should 
not expect A; and A,’ to be identical. Indeed, they cannot be. 
Either approach to the evaluation of the A’s and B’s is valid 
Thus it is not a mere coincidence that the final results obtained 
in Ref. 1 are in agreement with those of Ref. 2. 

3. The absence of a minus sign in the exponents of Eqs. (60) 
and (79)? is definitely an error of omission 

$. The final comment of Ref. (1) is concerned with the aero 
dynamic parameters, C, ™ and Ayw. It is stated! that Cr, the 
lift curve slope, is never negative, and that the damping param 
eter, A yw, is frequently negative. 

To begin with it is stated? that the graphical results which 
are presented are limited to the case of a conical body and the 


assumption of Newtonian aerodynamics. For these conditions 


Ci, 2(cos? B — sin? B) (a) 
1 1. \2 
c _ — 2cos? B (b) 
a 9 cos? B bret 
( 0 (c) A¢9 
(C ‘ r C, mm 7 
Amn a DI, Pp A ret Lr ef (d) 


Where B is the cone half-angle, J,.¢ is the cone length, /, is the 
static margin, and A,,; is the cone reference area 

\n examination of Eq. (7a) indicates that for B greater than 
15°, Cy, is negative, even for small angles of attack. Eq. (7b) 
indicates that for any value of B or 1,, Cy, is negative, and thus 
the damping parameter A yw is always positive. 

In general, the statement that CL, is never negative is incor- 
rect. For example, at angles of attack beyond the stall point, 
the lift curve slope of any aircraft is negative 

The original work on which Ref. 2 is based has recently been 
declassified.6 This report is reeommended for a more compre- 
hensive treatment of the problem of stability during re-entry 

Finally the writer would like to express his appreciation for the 
above comments. Their value is obvious in that points 1 and 3 
above were overlooked by both the author and the technical 
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A Note on the Compressible Turbulent 
Boundary Layer With Surface Mass Transfer* 


F. E. C. Culick 

Naval Supersonic Laboratory, 

Massachusetts Institute of Technology, Cambridge, Mass 
December 29, 1960 


Op POSSIBLE MEANS to alleviate the effects of aerodynamic 
heating is mass transfer or transpiration cooling. Practical 
considerations then require that one should investigate the be 
havior of turbulent boundary layers when gas is transferred 
through the surface. The most extensive analysis of the prob 
lem appears to be that of Rubesin and Pappas;' they fail, how 
ever, to find certain observed features of the effect of mass trans 
fer, although the influence of fairly high injection rates on heat 
transfer was described correctly Relatively simple calcula 
tions outlined here tend to agree with expected and observed 
results 

The boundary layer is regarded as comprising two crudely de 
fined regions: the “laminar sublayer’’ and an outer region in 
which turbulence is ‘‘fully developed.”’ Accordingly, two sets 
of equations are treated and the solutions matched at the inter 
face. Details of further assumptions are discussed in Ref. 2 
The results are reasonable and permit an approximate assessment 
of the effects of thermal diffusion 

Only the flat-plate problem is considered here. In the ab 
sence of detailed information a simple description of the flow in 
the turbulent region is obtained by invoking a form of ‘Rey 
nolds’ analogy,”’ assuming the ‘turbulent’ Prandtl and Schmidt 
numbers to be unity It follows as usual from the momentum 
and concentration equations for the turbulent region that the 
mean value of concentration, ¢, is linearly related to the velocity 
component uv. Then a “solution” to the energy equation is 
obtained by assuming the tempereture 7 to be a function ot wu 
Variations in the specific heat with ¢; are accounted for, and after 
satisfaction of the boundary corditions there results the useful 


relation: 


in which f = (Cys )/u,*; f df /di 
( 1+ arc — 1); u Us ue 


and hy is the normalized free-stream stagnation enthalpy C,,7 
u.2. Values at the edge of the laminar sublayer are denoted 
by ( ). Subscripts e and w denote evaluation in the free 
stream and at the wall 

The problem in the sublayer region is treated roughly as a 
perturbation of the case when the Prandtl number (Pr) and 
Schmidt number (Sc) are equal to one, justified in part because 
these quantities are always of order unity. Temperature dif 
ferences are assumed sufficiently small that Pr and Sc are func 
tions of c only. This condition is fulfilled by the available 


* This research was supported in whole or in part by the United States 
Air Force under Contract No. AF 49((38)-245, monitored by the Air Force 
Office of Scientific Research of the Air Research and Development Command 
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Fic. 1. Effeet of injection on wall concentration of helium 


measurements. It appears from results in Ref. 3 that a good 
approximation in the concentration equation is to neglect the 


effects of thermal diffusion so that this equation takes the more 


usual form 


pu(dc,/Ix) + pr(dc;/dy) 0/Oy[(u/Sc)(00,/0v)] (2) 


Then assuming ¢, a function of u yields a total differential equa 
tion for c;. This is further simplified by noting that in the sub- 
layer region the shear stress is approximately constant, especially 
at low injection rates. The resulting equation can be integrated 
formally and solved by a simple iteration procedure to give 


/Cw =1—0 = 1— (At — Bi?) (3) 
where A= J =. a= P : £= = wd 
1+, 1+, jr 2 
for which Sc varies linearly :° 
Sc = ji == Jo01 
In a helium-air mixture j; = 0.230 and j2 = 1.549. 


The energy equation can be handled in a similar manner, as- 
suming again that T is a function of uw. When these computa- 
tions are carried out, the Lewis number L = Pr/Sc is quite ac- 
curately represented by 
L = [ki/(ke — Kew)] + S; 

S = 0.210; K = (C,,. — 


k, = 0.852; kn = K — 0.295 
1 ery = 4.17¢1y 


Then the energy equation can be integrated formally, yielding 
integrals which can be evaluated numerically without difficulty. 

The influence of mass-transfer rate on skin friction and wall 
concentration of foreign gas may be calculated quite easily with 
a modified form of Covert’s work,! based on an approximate in- 
tegral method suggested by Liepmann.5 Details of these calcu- 
lations? will not be given here; the very simple formulas result: 


Cy he ( F\? F 


I; F Cy, ™ 
Clu : = F C ()) 
jz — Jj2 i] 
n+ (1+) Fic, +2 (14 : 2) 
21 2h Cy 
(pv)w 2 2F mr (PH )u 
~_ = and }\, = 
where (ea, C; C, Dick 


Cy, is the local skin friction coefficient with no mass transfer 
The result (4) agrees very well with the measurements of Ref. 6 
taken on a cone at a Reynolds number of 4 X 10% and M, = 3.24 
The data show virtually no changes in a Reynolds-number range 
2.93 X 108 to 5.80 X 108. However, a serious failing of Eq. 
(4) is that it does not predict the strong influence of Mach num 
ber observed in the reference. Values of wall concentration from 
Eq. (5) are shown in Fig. 1. 

Heat transfer (Stanton number) and the recovery factor are 
computed from the integrated energy equation and Eq. (1). It 
has been shown’ that with mass transfer at the surface and 
thermal diffusion included, the proper ‘‘adiabatic wall’’ bound- 


ary condition is 
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[x(OT'/Oy)] au [pou RTRT/c1))a 


which can be written in terms of previously defined quantities as 


lea = —Ghas; —G 21, [(RPrkt)/(oCp)] 6) 


where kt is the thermal-diffusion ratio and R the gas constant 
for the mixture. The recovery factor and Stanton number are 


then given by 


1 , GW 
} : rF+2Z 
(1 — GW) (y — 1M? 


1— GW 
1+ GW/(r'/2)(y — 1)M?2 + [1 — (Tx/T-)] 


with 


W | dé me cil — =A : l ke — Ke) 
’ x Xl P) 

and r’, St’ are the values with mass transfer but no thermal dif 

fusion.? For simplicity, effects of the latter appear here only 

through the boundary condition (6). 

Fig. 2 shows some computed results for a helium-air mixture, 
and also the special case of air injection. The velocity at the 
sublayer edge has been left unknown so far as its variation with 
injection rate is concerned. Evidently the calculations are not 
excessively sensitive to its value. The dashed curves including 
the ‘‘correction’”’ for thermal diffusion cannot be regarded as 
especially accurate, but the obvious trend appears to be true 
The overshoot in 7/ro at low injection rates in Fig. 2 can be 
ascribed to the factor kt/c, in C; slight changes in this ratio can 
change the unsatisfactory behavior. There appears to be some 
uncertainty in the values of kt. 

For the most part, the results are satisfactory, with two ex- 
ceptions: the calculated skin friction does not exhibit the ob- 
served dependence on Mach number, and the predicted recovery 
factor lies considerably below measured values. It is perhaps 
an important result, although not truly conclusive, that if thermal 
diffusion is included, then the latter error is reduced while the 
heat transfer is only slightly affected. This might have been 
anticipated from estimates in Ref. 3 (pp. 76-80). 
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Fic. 2. Comparison of theory and experimental results for re- 


covery factor and Stanton number: helium injection. 
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The Epicyclic Analysis as Applied 
to Free-Flight Data 


Peter J. Mantle 

Graduate Student, Guggenheim Aeronautical Laboratory 
California Institute of Technology, Pasadena, Calif. 
December 27, 1960 


SYMBOLS 


K;(j = 1,2,3) = rotating vectorsin amplitude plane 
k = pitch radius of gyration 
ly = reference length 
m = body mass 
p = steady rate of roll 
Sr = reference area 
l = forward velocity of body c.g 
Iz = rolling moment of inertia about body axis 
I = pitching moment of inertia about body axis 
Ax = horizontal displacement between peaks in amplitude 
plane 
a = pitch incidence relative to rolling body axes 
a@ = pitch incidence relative to nonrolling coordinates 
8 = sideslip angle relative to rolling body axes 
B = sideslip angle referred to nonrolling coordinates 
u = complex angle of attack (8 + ia) 
wu} = ‘‘total’’ angle of attack V a? + Bp? 
j(9 = 1, 2) = damping parameters in epicyclic solution 
w3(j = 1, 2) = angular rates of vectors K; 


Subscript 

1 = evaluated at in-phase position —i.e., at peaks 

rw A COMPLEMENTARY NOTE to Ref. 1 a short discussion is 

presented on the application of the epicyclic analysis to 
free-flight data of axisymmetric bodies, with particular reference 
to data obtained in an aeroballistic range. 

After Fowler et al.2 demonstrated the epicyclic behavior of 
spinning shells, Nicolaides* introduced the technique of analyzing 
the data in the complex amplitude plane where the coordinate 
and abscissa are pitch incidence (&) and sideslip angle (8), 
respectively. Nicolaides’ analysis was characterized by the use 
of a nonrolling body-axis system which resulted in the complex 
angle of attack being represented by the sum of three rotating 
vectors (the tricyclic form). If, however, by a rotation of axes 
we return to rolling coordinates, where the rotation is given by 
the integral of the roll rate, then the epicyclic form of the solution 
is obtained and we may write* 


B+ ia = Kyl + or — po) 4 KX, ele + Hor — pot 4 K, (1) 


The quantities K,, Ke, K; are the initial values of the complex 
angle of attack, and specifically, Ks; is the trim vector associated 
with body asymmetries and control inputs. The resonance 
phenomenon related to this trim vector has been considered by 
Murphy.® 

The epicyclic analysis as presented here results in a particularly 
simple data-reduction technique and avoids the iterative pro- 
cedures required in both Nicolaides’ tricyclic analysis and 
Danis’ epicycle method. The method is best described with 
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the aid of anexample. For this purpose the flight-data obtained 
from the trajectory of a flared-base projectile in the aeroballistics 
range at the Canadian Armament Research and Development 
Establishment (CARDE) will be used. The flight data of a 
typical trajectory are presented in Fig. 1 

Each experimental point in this amplitude plane corresponds 
to the value of the total angle of attack measured at a particular 
station. Proceeding from point to point along the curve describes 
the flight history—i.e., Fig. 1 represents the trace of the nose of 
the projectile during its flight downrange. The value of Y at 
each point has been omitted for clarity. The development of 
the method is given in Ref. 4 but the pertinent results are sum 
marized here. The technique is based on the condition that the 
loops in the amplitude plane (Fig. 1) correspond to the in-phase 
positions of the vectors K, and Ky, (note the vector K; can be 
eliminated by a simple translation of the origin). From this the 
following important results are obtained which provide solutions 
for the static stability and damping parameters 


Cm — [(w/Ax)? — (Iz/I)*%( po/ U)?|(21/pSA,) 2) 
Cya — (lr R)*( Cm, + Cm) —(2m/pS,)(v + v2); (3) 
In Eq. (3) the result is limited by the condition that << 1 


which is satisfied by most conventional configurations. The 
value of the steady roll rate can be obtained from the rate of 
migration of the loops in the amplitude plane, but it is more 
usually obtained from the original free-flight data. For slender 
bodies rolling at moderate rates the second term in Eq. (2) 
can be neglected, giving the simple approximate value for the 
static stability coefficient 

Cm = —(n/Ax)%(2I/pS,l,) (4) 


The reduction technique utilizes two quantities taken from the 
amplitude plane (Fig. 1). These are X,, the distance flown to 
the peaks of each loop, and |u!», the total angle of attack meas- 
ured ateach peak. Plotting Y, versus the peak number (counted 
from the first convenient peak in the trajectory) gives a straight 
line of slope AY which is shown, for the cited example, in Fig. 2. 
The slope of this line and Eq. (2) or (4) give the value of the 
static-stability parameter. 

Plotting the value of |u| at each peak on a logarithmic scale 


against the linear distance to peak gives Fig. 3. The slope of 
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Fic. 1. Amplitude plane, rolling coordinates. 
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Fic. 2. (Left) Determination of static stability (C,,,). Fie. 3. 
(Right) Determination of total damping. 


this straight line (1) + v2) is directly proportional to the damping 
of the system as given by Eq. (3). This completes the reduction 
technique. 

This epicyclic method has been applied successfully to many 
free-flight data obtained on axisymmetric bodies in the aero- 
ballistic range at CARDE. These data pertained to configura- 
tions tested over a large range of flight Mach numbers and roll 
rates. Although this method is only applicable to the case of 
steady roll, it is obvious that the trajectory can be divided into 
regions over which a mean value of the roll rate is taken and the 
solutions superposed. 

In conclusion, it is seen that this method is particularly simple 
to apply and avoids the iterative procedures previously em- 
ployed and also the need for digital computers. A restriction 
on the method is that only a small portion of the available ex- 
perimental data is actually used, which places a more stringent 


requirement on the accuracy of such data. 
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Shock Propagation in Time Department Ducts* 


Roy Gundersen 

Associate Professor of Mathematics, Illinois Institute of Technology, 
Chicago, III 

December 29, 1960 


LTHOUGH problems with ducts of time-dependent cross- 
sectional area are rare (one example would be a boundary 
layer growing with time), it is of interest that the pressure change 
behind a shock propagating through such a duct is given by the 
solution for the time-independent duct by making a simple 
change of the time variable. 
The equations which govern the one-dimensional, unsteady 
flow of an inviscid, ideal compressible gas in a tube of varying 


cross-sectional area are: 


9 > 
2 e* 
Uy + UUz + CC; = St 
7-1 Vy — 1)e, 
» ») 
Ec; + Euc, + cEu, + cE; + ucE, = 0 
y¥—-1 y- 1 
Sit usz = O 


* This work was supported by a grant from the National Science Founda- 


tion 
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where u, c, s, ard E are, respectively, the particle velocity, the 
local speed of sound, the specific entropy, and the cross-sectional 
area. Partial derivatives are denoted by subscripts and  ¥ is 
the ratio of specific heats at constant pressure c, and at constant 
volume c, 

If the cross-sectional area of the tube is expressed in the form 
E(x,t) = Ey + E,(x,t), where E, is the original uniform area and 
E\(x,t) is the perturbation, and the above equations are linearized 
in the neighborhood of a known (isentropic) solution denoted by 
the subscript zero, the following system of linear equations, which 
govern the first approximation to the behavior of the disturbed 
flow (not assumed isentropic) is obtained: 


i ++ 1 1-3 
A, + (tm + o)A; + = { + = B la; 


—(1/2E5)( toc Ei, — < E},) + (1/2) H 


l 
3-34 y+ 1 
Bi + (tt — c)Br + a 
—(1/2Fo)(tcoFir, + coF\,) — (1/2)H | 
Uo Co uo . a 
where 9 + —e a, _ 2 + eT = 8 
1; Cc} ch Cc} 
4 1, al B 
2 y¥- 1 2 y¥- 1 


and Hocry(y — 1) = co? 
The derivation closely parallels that of Ref. 1, which may be 


consulted for details 
For the case of an initially uniform flow, a and 6 are constants 


and the general solution to the system (1) is: 


Hy = 2R'(x — wuol) 


A = Flx — (to + co)t] + co R(X — wt) — 
E,[x, (to + co)t/uo] uoco/2Eo( to + ci b (2) 
B = G(x — (uo — colt] + co R(X — wt) — 
E, (x, (to — €o)t/tto] Uoco/2Eo( tt — co) 


When an originally uniform shock wave propagates through a 
tube of varying cross-sectional area, the shock is perturbed, the 
shock strength is altered, and the subsequent flow is nonisen 
tropic. There are two contributions to the disturbance—viz., 
a permanent perturbation due to the area variations and a 
transient perturbation propagated downstream with sonic ve- 
locity relative to the flow behind the shock, due to reflections 
from the shock of the permanent perturbation 

Suppose the area variations are restricted to the region x > 0 
and that the shock of arbitrary strength initially moves with 
constant velocity through the region x < 0. If the shock is 
generated by a piston moving with constant velocity, it is as 
sumed that the piston is always sufficiently far to the left so that 
no reflections come back from it to interfere with the interaction 
under consideration. 

Let the flow in front of the shock be denoted by the subscript 
1 and the flow behind by the subscript 2. Let P denote the pres 
sure, w the shock velocity, v w— &, Mm u/c, M = v/c so 
that the following relations hold: 


U2 = uw + 2Mic(1 — W,7?)/y + 1 

Ci / C2 = (yr + 1) Ms Mi(4 — 1+ 2),~?) 

P2/P, = 2yM,2/(y + 1) — (y — 10/9 +-:1) 

M2 = (7 — 1+ 2M,~?)/[2y — (vy — 1).7?3] (3) 
my» = (y + 1)Mom/Mi(y — 1 + 2472) + 


2(1 — M2?)Moxy + 1) 
Mie(y — 1 + 2M,~*) 


(y + 1) 


) ») 


The present problem is identical with that solved in Ref. 2 
save for the added dependence of E; on ¢t. Hence the results 
may be given with only minor modifications. The perturbation 
of the shock will be obtained by the use of the general solution, 
Eqs. (2), by setting F = 0, since, from the formulation of the 
problem, the contribution given by F is not included, for there 
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is no mechanism downstream of the shock which could give rise 


—(y + 1) P2 — Pi)th/Mia(l — M72) 1 + 2M. + M72) 4 
+ 1) Ps — Pi) + M,~?)P,/2yP,(1 — Me? 1 + 2M. + M, 
+ 1)m(P2 — P,)E;[x, (1 + me)t/me) /Ep 1 + me)My(1 -— Mi, 


> 


ro complete the present solution, the arbitrary function G of Eq 
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to such a term. Thus, denoting perturbations of a base quan 


tity by a bar, setting F = 0 in the second of Eqs. (2), and utilizing 


the perturbation form of the relations (3) we obtain 


2) — 2 P2 — P;)&/q(1 ML + M - 

1+ 2M.+ M,~?) - 
P, — P,)E, |x, (1 + me)t/m Ewl + n 1+2M.+ MV }) 
2) must be determined. This is most easily obtained by adding the 


last two of Eqs. (2), evaluating the result on the shock, x wt, and replacing P» by its equivalent from Eq. (4). This yields 

rn) /e —ti|[rrA, A/c 1 + AMe)| 201 — Me?)/Miq(1 — 14,7-2)(1 + 2Me + M72) - 

1i[7A, A/c 1 + Me)|(1 — Me?)/(y + 1)g(1 + 2Me + M2) 4 

Pi[rr, X/cA1 + M2)](1 + Mi-2)/yPi(1 + 2Me + M,-2) + meE,[rd, (me — 1)d/moce( Me + 1))/2E( me — 1) 4 

) m2/2(1 + me) — 2my(1 — Me?)/(1 + me)M\(1 — MW, 7-?)1 + 2Me + 7?) — 401 -— Me? y¥ + 1)(1 + me)1 + 242. + M ‘ 

) Ey [| rd, (mm L)A/mtocol Mo + 1 Ey; 5) 

where 7 Wo + m 1+ Jz) and Av x — (tt — c2)f 

rhus, the complete solution for Pris 
Ps/(P2 — P;) —)F,lx, (1 + mie)t/me| /(me + 1) + Eylx, (me — 1)t/me) /(me -— 1) - 

E,|rv\, (2 — 1)XN/nt0¢0 Meo + 1 ms — 1)4(y + 1)me/4E (1 — Me?) + 

(Cy + L)mte/4(1 + me)(1 — Me?) — (Cy + 1) /(1 + me)M\(1 — M,7?)(1 + 212 + 14,7?) 4 


—2/(1 + mel + 2Meo + W,7*){ Ei rd, (me + 1)A/miece Me + 1 
@:2/c(1 + 2Mo + M,~2) Py + 1)1 + -?)/2yPi(1 + 2M 4 


where the omitted arguments of @, ¢; and P, are the same as in Eq. (5 
In the region directly behind the shock, the pressure change is 
given by Eq. (4). When the flow behind the shock is subsonic, 
i disturbance will be reflected downstream of the shock, and the 
When the flow behind 


the shock is supersonic, this disturbarce will propagate upstream 


pressure perturbation is given by Eq. (6) 


of the shock, and the pressure perturbation is given by the entire 
Eq. (6) 
It is easily seen that when there is no time dependence ir the 


irea variations, the present results reduce to those presented in 
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Stagnation-Point Heat Transfer in Equilibrium 
Dissociated Air With Variable Prandtl 
and Lewis Numbers 


Aemer D. Anderson* 
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Sunnyvale, Calif 

January 3, 1961 


SYMBOLS 
= mass fraction of component 1 
= specific heat at constant pressure of component i 
Cp = Ci Cpi 
t 
= total specific heat at constant pressure of mixture 


D = diffusion coefficient 
hi® = heat of formation of component i at absolute zero 
h = total enthalpy of mixture } f 4 cpi dT + hi 
é 
h = chemical enthalpy at boundary layer edge 2 Cie hi 
2 
k = ordinary thermal conductivity due to molecular collision 
ke = total equilibrium thermal conductivity 
Le = Lewis number, pDCp/k 
p = pressure 
Pr = Prandtl number ulp/k 
2 = heat flux 
T = absolute temperature 
u = velocity in x direction 


* Thermodynamicist 


M,-?)(1 — M27? 6) 


velocity in y direction 

x surface coordinate 
normal coordinate 

p = density 


Mu absolute viscosity 


Subscripts 


e = boundary-layer edge 


- = wall 


ern THEORY of stagnation-point heat transfer in dissociated 
air was set forth by Fay and Riddell in 1958.!. They re 
ported numerical solutions to the stagnation-region boundary 
layer equations for air in equilibrium, which were obtained as 
suming constant Prandtl and Lewis numbers in each solution 
These calculations were for P? 0.71 and 10 < Le < 2.0 
Their correlation of these solutions may be written as follows: 


g = 0.76 Pr~-* ( pee du, /dx) ( prptte/ pete)! X 


[1 + (Le®%5? — 1)ha/h he — In) (1) 
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TABLE 1 





Numerical Solutions 

















t | | | = | 
de % a nwhing| Dtu/lo pont a sa = | ten 
— + — — <3 + + 
0.1 3500 500 2935 | a7 2.56 | 0.77 1.32 | 0.34 | 1.00 
0.1 6000 500 | IH 217 | 3.54 0.73 1u | 0.54 | 0.98 
1.9 2500 500 | 1261 | 217 | 1.92 | 0.76 | 1.0 | 0,02 | 1,00 
1.0 2500 1000 1261 Lso 1.49 0.77 1.40 0,02 1,00 
1.0 3500 | 500 | 2477 | 217 | 2.43 | 0.76 1.% | 0.2u | 0.99 
1.0 | 3500 | 1000 | 2477 | 450 | 1.89 | 0.77 | 1.¥ | 0.2L | 1.01 
1.0 } 5500 | 500 | 5140 | 217 | 325 | 0.75 1.26 | 0.39 | 0.98 
1.0 | $500 1000 | SIMO | U50 =| 2653 | 0.75 | 1.26 0.39 | 1,00 
1.0 | 5500 2000 | 5140 | 967 | 1.69 | 0.75 | 1.26 , 0-39 | 1.03 
1.0 | 5500 | 3000 51 1705 1.53 0.75 1.26 39 1.03 
1.0 | ss00 | wooo | Simo | 334 | 1.20 | 0.73 | 1.17 | 0.39 1.04 
1.0 6500 500 8511 217 3.78 0.73 1. 1.56 97 
1.0 6500 000 8511 450 2.9% 0.73 1.14 0.5 1,90 
1.0 70K 500 11,140, 217 4.12 0.73 1.07 64 0.98 
1.0 | 7000 1000 | 11,1,0) 50 | 3.20 | 0.73 | 1.07 | 06h | 1.00 
1.0 7900 200 1,140 967 2.37 0.73 1.07 0.41 1.03 

00 3000 11,140; 1705 1.93 0.73 1.07 0.54 1.04 
i. 7900 4000 11,10, 331k 1,51 0.72 0.9 0.64 1.07 
1.0 8000 500 15,390 217 4.68 0.73 1,00 0.59 0.97 
1.0 3009 1000 16,390, 450 0.73 1.90 0.69 1.00 
10.9 4000 500 2737 217 217 2.40 421 0,99 
13.0 700 500 7552 217 3.73 0.74 1.21 O47 0.% 
100,90 Sao 500 2550 217 20% 0.77 1.43 0.16 0.99 
190.0 8000 500 80L9 17 3.90 0.74 1.22 oP % 


This expression defines the effect of the level of the density- 
viscosity product and the effect of its distribution through the 
boundary layer, at least for the Sutherland viscosity model used. 
It also indicates the influence of the level of the Prandtl and 
Lewis numbers. 

In dissociated-air boundary layers the Prandtl number 1s in 
fact nearly constant; this is not true, however, of the Lewis 
number.2* The effect of Lewis-number distribution through 
the boundary layer is left undefined by Eq. (1) and could very 
well be important at high enthalpy levels. The purpose of this 
note is to show that Eq. (1), with its attendant simplicity, can 
be used in situations where the Prandtl and Lewis numbers do 
vary across the boundary layer, when representative values of 
these parameters are selected in a particular way. 

The present uncertainty in the determination of the transport 
properties of dissociated air is considerable.” The calculations 
of Hansen? were used in this study since they provide, in con- 
venient form, a complete model for both the thermodynamic and 
the transport properties. 

The solutions reported herein were carried out, primarily, for 
a pressure of l atm. The viscosity given in Ref. 2 for this pres- 
sure is shown normalized with the Sutherland viscosity in Fig. 1. 
The Prandtl and Lewis numbers from Ref. 2 are also shown in 
Fig. 1 for this pressure. These property values are not strongly 
pressure-dependent in the enthalpy and pressure ranges of in- 
terest here. 

The total equilibrium thermal conductivity, k;, as tabulated 
by Hansen, describes the overall energy transfer in the presence 
of a temperature gradient. This includes the energy transfer by 
diffusion as well as that by ordinary conduction due to molecular 
collision. The total specific heat is defined as 

Coe = (Oh/OT)p 
where /: includes both the sensible and the chemical enthalpies. 

Using these properties, the boundary-layer conservation equa- 

tions may be written for the stagnation region, as follows: 


Mass: 
(0/O0x)(pux) + (0/dy)(pux) = 0 (2) 
Momentum: 
Ou Ou re) Ou du, , 
pu + pv = + Pelle (3) 
Ox ‘ oy oy 4s oy dx 
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Energy: 
oh oh O (fk ch 


pu 5 01 = +) 
x 


s Oy OV \Cp: OY 
Numerical solutions of these equations were obtained using the 
transformations and calculation procedure described in Ref. 4 
A revision was made so that the effect of the compressibility, Z, 
was included in the calculation of the density. 

The boundary conditions for which the various solutions were 
carried out are shown in Table 1. Also shown are the param 
eters used in the Fay-Riddell correlation [Eq. (1)| and a com 
parison of the correlation with the numerical solutions 

The quantity Le is the average of the maximum and minimum 
values of the Lewis number in the boundary layer. The repre 
sentative Prandtl number, Pr, was similarly evaluated With 
this choice of these parameters, Eq. (1) is seen to correlate the 
numerical solutions within +5 percent, with one exception at a 
very high wall temperature. 


This agreement is somewhat to be expected since the use of 


either the maximum or minimum Lewis number in Eq. (1) 
would change the prediction by only 15 percent at the most 
(increasing the scatter to +20 percent). These solutions do 
show, however, that neither the Lewis number nearest the wall 
nor the Lewis number at the edge of the boundary layer (where 
the chemical enthalpy is largest) is more representative of the 


entire boundary layer. 
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Comparison of Flow-Direction Probes 
at Supersonic Speeds 
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eS working on flight instrumentation, wind 
tunnel instrumentation and aircraft air-induction-system 
controls are concerned with measuring Mach number and flow 
direction. The selection of probes for measuring Mach number 
has been discussed.! Various wedge, conical, hemispherical and 
pitot probes which provide a pressure signal related to flow 
direction are compared here on the basis of theoretical and ex- 
perimental sensitivity (rate of change of pressure difference or 
of pressure ratio with flow-direction angle) data. 


THEORY 
Static-pressure taps on a conical probe have been used to 
measure flow direction.2. The theoretical sensitivity, if the flow 
is aligned with the cone axis, is 
AP/Qe = 0.034906(P/Q)n/P (1) 
where AP is the difference between pressures at two diametrically 
opposite taps, Q is the free-stream dynamic pressure (7~Ps0M?/2), 
e is the flow angle in the plane of the taps in degrees, and n P 
and P/Q may be computed from tables.* The sensed pressures 
are larger for a probe with pitot tubes near the cone surface 
instead of static taps, and the sensitivity is 


AP/Qe = (0.034906P,/Q)(n/P + KX) (2) 
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Fic. 1. Sensitivity of probes with static taps. 


where the pitot pressure P, is computed from P and M,3, X is 
tabulated,‘ and 
: 8 (2M,2 — 1)y}1 + [(y — 1)/2]M,7}?? 
K (3) 
Ws «4 M(2yM,Z — y + 1) 
if the cone surface Mach number /, exceeds one. 

The theoretical sensitivity of a two-dimensional wedge probe 
at zero flow angle may be derived from an extension of the anal- 
ysis of Ref. 1. In terms of the nomenclature of Ref. 1, the 
sensitivities are, for two surface static taps, 

AP/Qe = —(4t/yM*c)ca/c (4) 
0.069813 Vie — ] 
for a flat plate (w = 0); and, for two pitot tubes between the 
oblique shock and the surface, 
AP Ve — (4 4 M* )[d( ca = diniMa (5) 
by{ 1+ - M? 
0.069813 } 2 
V M2 —1 yM 
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COMPARISON OF PROBES 


The theoretical sensitivities of wedge and conical probes using 
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Fic. 2. Sensitivity of probes with pitot tubes 
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Fic. 3 Effect of geometry on probe sensitivity 


surface static taps and an empirical correlation (Ref. 5) for 
hemispherical probes are shown in Fig. 1. The theoretical sensi 
tivities of probes using pitot tubes are shown in Fig. 2. Meas 
urements with pitot tubes on a 26° included-angle cone for 
Hamilton Standard Division showed slightly lower sensitivities 
than were computed from Eq. (2). Comparison of Fig. 2 with 
Fig. 1 shows that greater sensitivities can be obtained with pitot 
tubes than with static taps on probes at Mach numbers exceeding 
roughly two. The measured sensitivities of a hemisphere with 
two short, projecting, pitot tubes at 45° to the axis (similar to 
the Douglas-Giannini probe®) and of a Y-shaped probe with an 
included angle of 120° between the two pitot arms are also shown 
in Fig. 2. The sensitivities of the latter probe were derived 
from tests of a pitot tube in the UAC Research Laboratory super 
sonic wind tunnel at flow angles from 5° to 80 

The variation of the sensitivity of several probes with included 
angle is shown in Fig. 3 for a Mach number of 3.0. The sensi 
tivities of pitot tubes on a probe are a maximum for a 10° wedge 
and for a 28° cone and become negative at large included angles 
The sensitivities of static taps on wedges and cones increase with 
included angle. The sensitivity of two pitot tubes is largest for 
an included angle of roughly 110° between the arms 

The probes discussed above are compared in Fig. 4 for appli 
cations where a sensed pressure ratio rather than a pressure 
difference is used to determine flow direction. As shown by the 
line labeled ‘‘typical locus of constant accuracy,” friction and 
threshold effects in the sensor require that a probe providing 
signals with a lower pressure level have a larger sensitivity in 
order to maintain constant accuracy in sensed flow direction 
Therefore, static taps on a 60° wedge can provide the same 
accuracy as pitot tubes on a flat plate 

Flow-direction angle in one plane can be determined at super 
sonic speeds by measuring a pair of pressures from a flow-direc 
tion probe. Selection of an appropriate probe geometry will 
optimize the pressure signal and thereby minimize errors in the 


derived flow direction 
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Fic. 4. Effect of geometry on pressure-ratio sensitivity and 
pressure level. Number beside symbol is included angle of 
probe (2w) in degrees. 
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action with the surrounding medium (air, fluid); this 
effect will be discussed in a future publication. Also, 
another limitation of the present analysis should be 
indicated: since the strains and stresses in the z- 
direction have been neglected, the results do not seem 
to be applicable for frequencies and modes for which the 
distance between nodal lines is smaller than about twice 


the thickness of the plate. 
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THE CHASE MANHATTAN BANK 

CHICAGO AERIAL INDUSTRIES, INC. 

CONTINENTAL MOTORS CORPORATION 

CORNELL AERONAUTICAL LABORATORY INC. 

COX & COMPANY, INC. 

CURTISS-WRIGHT CORPORATION 

DANIEL, MANN, JOHNSON, & MENDENHALL 

THE DECKER CORPORATION 


DEFENSE ELECTRONIC PRODUCTS DIVISION, RADIO 
CORPORATION OF AMERICA 


DEFENSE SYSTEMS DIVISION, GENERAL MOTORS 
CORPORATION 





DEL MAR ENGINEERING LABORATORIES 

DELCO-RADIO DIVISION, GENERAL MOTORS CORPO- 
RATION 

DELCO-REMY DIVISION, GENERAL MOTORS CORPO- 
RATION 

DOAK AIRCRAFT COMPANY, INC. 

DOUGLAS AIRCRAFT COMPANY, INC. 

DZUS FASTENER COMPANY, INC. 

EASTERN AIR LINES, INC. 

EATON MANUFACTURING COMPANY 

EDO CORPORATION 

ELASTIC STOP NUT CORPORATION OF AMERICA 

ENGINEERING SUPERVISION COMPANY 


FAIRCHILD CAMERA AND INSTRUMENT CORPORA- 
TION 


FAIRCHILD STRATOS CORPORATION 

THE GARRETT CORPORATION 

GENERAL APPLIED SCIENCE LABORATORIES, INC. 
GENERAL DYNAMICS CORPORATION 

GENERAL ELECTRIC COMPANY 

GENERAL PRECISION, INC 

THE B. F. GOODRICH COMPANY 

GOODYEAR AIRCRAFT CORPORATION 
GRUMMAN AIRCRAFT ENGINEERING CORPORATION 
HOFFMANN AIRCRAFT COMPANY 

HUGHES AIRCRAFT COMPANY 


INSURANCE COMPANY OF NORTH AMERICA COM- 
PANIES 


IBM SPACE GUIDANCE CENTER 

THE INTERNATIONAL NICKEL COMPANY, INC. 

ITT FEDERAL LABORATORIES, A DIVISION OF INTER- 
NATIONAL TELEPHONE AND TELEGRAPH COR- 
PORATION 


ANITROL eae DIVISION, MIDLAND-ROSS COR- 
PORA 


JOHNS-MANVILLE SALES CORPORATION 
KOLLSMAN INSTRUMENT-CORPORATION 
LAVELLE AIRCRAPT CORPORATION 
LEAR, INCORPORATED 

C. W. LEMMERMAN, INC. 

THE LIQUIDOMETER CORPORATION 


ARTHUR D. LITTLE, INC. 








LOCKHEED AIRCRAFT CORPORATION 
THE MARQUARDT CORPORATION 

THE MARTIN COMPANY 

McDONNELL AIRCRAFT CORPORATION 
MELETRON CORPORATION 
MINNEAPOLIS-HONEYWELL REGULATOR COMPANY 
MOBIL OIL COMPANY 

NATIONAL CREDIT OFFICE, !NC 

NORTH AMERICAN AVIATION, INC. 
NORTHROP CORPORATION 

PAN AMERICAN WORLD AIRWAYS, INC. 
THE RALPH M. PARSONS COMPANY 
PIASECKI AIRCRAFT CORPORATION 
POLY INDUSTRIES, INC 

REPUBLIC AVIATION CORPORATION 
ROHR AIRCRAFT CORPORATION 

RYAN AERONAUTICAL COMPANY 
SHELL OIL COMPANY 

SIMMONDS PRECISION PRODUCTS, INC. 
SOLAR AIRCRAFT COMPANY 
SOUTHWEST PRODUCTS COMPANY 
SPACE TECHNOLOGY LABORATORIES, INC. 
R. DIXON SPEAS ASSOCIATES 


SPERRY GYROSCOPE COMPANY, DIVISION OF SPERRY 
RAND CORPORATION 


STANDARD-THOMSON CORPORATION 
STANLEY AVIATION CORPORATION 
THERM, INCORPORATED 

THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC. 
TRANS WORLD AIRLINES, INC. 


UNION CARBIDE METALS COMPANY, DIVISION OF 
UNION CARBIDE CORPORATION 


UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 

UNITED STATES AVIATION UNDERWRITERS, INC. 
VEHICLE RESEARCH CORPORATION 

WESTERN GEAR CORPORATION 


WESTINGHOUSE ELECTRIC CORPORATION 




















